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Abstract. The Numerical INJection Analysis (NINJA) project is a collabora- 
tive effort between members of the numerical relativity and gravitational-wave 
data analysis communities. The purpose of NINJA is to study the sensitivity of 
existing gravitational-wave search algorithms using numerically generated wave- 
forms and to foster closer collaboration between the numerical relativity and data 
analysis communities. We describe the results of the first NINJA analysis which 
focused on gravitational waveforms from binary black hole coalescence. Ten nu- 
merical relativity groups contributed numerical data which were used to generate a 
set of gravitational-wave signals. These signals were injected into a simulated data 
set, designed to mimic the response of the Initial LIGO and Virgo gravitational- 
wave detectors. Nine groups analysed this data using search and parameter- 
estimation pipelines. Matched filter algorithms, un-modcllcd-burst searches and 
Baycsian parameter-estimation and model-selection algorithms were applied to 
the data. We report the efficiency of these search methods in detecting the nu- 
merical waveforms and measuring their parameters. We describe preliminary 
comparisons between the different search methods and suggest improvements for 
future NINJA analyses. 
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1. Introduction 

Binary systems of compact objects, i.e., black holes and neutron stars, are among the 
most important objects for testing general relativity and studying its astrophysical 
implications [1]. The general solution of the binary problem in Newtonian gravity 
is given by the Keplerian orbits. In general relativity, the Keplerian orbits for a 
bound system decay due to the emission of gravitational radiation, leading eventually 
to the merger of the two compact objects and to a single final remnant [2-4]. 
The decay of the orbits is due to the emission of gravitational waves and these 
waves carry important information about the dynamics of the binary system. In 
particular, the waves produced during the merger phase contain important non- 
perturbative general rclativistic effects potentially observable by gravitational-wave 
detectors. Gravitational waves could be detectable by the current generation of 
gravitational wave detectors such as LIGO and Virgo [5,6], and detection is very 
likely with future generations of these detectors. 

Two important advances have occurred in recent years that have brought us 
closer to the goal of observing and interpreting gravitational waves from coalescing 
compact objects. The first is the successful construction and operation of a world- 
wide network of large interferometric gravitational-wave detectors; these include the 
three LIGO detectors in the United States, Virgo in Italy, TAMA in Japan [7] and the 
GEO600 detector in Germany [8]. The TAMA detector was the first interferometric 
detector to achieve its design goals, and it collected science data between 1999 and 
2003 [7]. The LIGO detectors started observations in 2002 [9]. From 2005 to 2007 
these detectors operated at design sensitivity collecting more than a year of coincident 
data from the three LIGO detectors; these observations are referred to as the "fifth 
science run" (S5) [10]. The Virgo detector is also close to achieving its design goals 
and collected six months of data coincident with the last six months of the LIGO 
S5 run (referred to as VSR1) [11]. The GEO600 detector has been operating since 
2002 in coincidence with the LIGO instruments [8] . The two 4km LIGO detectors are 
currently being upgraded to improve their sensitivity by a factor of 2-3 (Enhanced 
LIGO [205]) and will resume observations in 2009. Upgrades to the Virgo detectors to 
yield comparable sensitivity to Enhanced LIGO arc proceeding on a similar schedule. 
During this time, the GEO600 and the LIGO Hanford 2km detector continue to make 
best-effort observations (called "astro-watch" ) to capture any possible strong events, 
such as a galactic supernova. Following the Enhanced LIGO and Virgo observations, 
the Advanced LIGO [206] and Virgo [207] upgrades will improve detector sensitivities 
by a factor of ~ 10 above the Initial LIGO detectors; these upgrades are expected 
to be complete by 2014. There are also plans to build a second-generation cryogenic 
detector in Japan known as LCGT [12]. Searching data from these detectors for 
weak gravitational wave signals over a vast parameter space is a challenging task. 
The gravitational- wave community has invested significant resources in this effort. A 
number of searches on S5/VSR1 data for un- modelled bursts and binary coalescence 
are in progress and many results, including those from previous science runs, have 
already been reported [13-23]. 

The second important advance has been the impressive success of numerical 
relativity in simulating the merger phase of binary black hole (BBH) coalescence. 
The first breakthroughs occurred in 2005 with simulations by Pretorius [24], closely 
followed by the independent Goddard and Brownsville (now at RIT) results [25,26]. 
Since then, a number of numerical relativity groups around the world have successfully 
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evolved various configurations starting from the inspiral phase all the way through 
the merger to the final remnant black hole (for recent overviews on the field see 
e.g. [27-29]). This has led to important new physical insights in BBH mergers. These 
include the prediction of large recoil velocities produced by asymmetric emission 
of gravitational radiation during the merger process [30-48] and the prediction of 
the parameters of the remnant Kerr hole for a wide class of initial states [32,49- 
61]. Since the inspiral, merger and coalescence of black holes are also among the 
most important targets of gravitational-wave detectors, we expect that the detailed 
information provided by numerical simulations can be used to increase the reach and 
to quantify the efficacy of data analysis pipelines. Indeed the driving motivation of 
research on numerical simulations of black-hole binaries over the last few decades has 
been their use in gravitational- wave observations. 

Thus far, most searches for gravitational waves from BBH mergers have relied on 
post-Newtonian results, which are valid when the black holes are sufficiently far apart. 
Within its range of validity, post-Newtonian theory provides a convenient analytic 
description of the expected signals produced by binary systems. The numerical 
relativity results, on the other hand, have not yet been synthesised into an analytic 
model for the merger phase covering a broad range of parameters, i.e., a wide range 
of mass ratios, spins and if necessary, eccentricity; there has however been significant 
progress for the non-spinning case [51,62-72]. Similarly, despite significant progress, 
there is not yet a complete detailed description over the full parameter space of how 
post-Newtonian and numerical simulations are to be matched with each other. On 
the data analysis side, many pipelines, especially ones that rely on a detailed model 
for the signal waveform, have made a number of choices based on post-Newtonian 
results, and it is important to verify that these choices are sufficiently robust. More 
generally, it is necessary to quantify the performance of these data analysis pipelines 
for both detection and parameter estimation. This is critical for setting astrophysical 
upper limits in case no detection has been made, for following up interesting detection 
candidates, and of course for interpreting direct detections. Work on this to date 
has primarily used post-Newtonian waveforms. Numerical relativity now provides an 
important avenue for extending this to the merger phase. 

There are significant challenges to be overcome before numerical relativity results 
can be fully exploited in data-analysis pipelines. The Numerical INJection Analysis 
(NINJA) project was started in the spring of 2008 with the aim of addressing these 
challenges and fostering close collaboration between numerical relativists and data 
analysts. Participation in NINJA is open to all scientists interested in numerical 
simulations and gravitational-wave data analysis. NINJA is the first project of its 
kind that attempts to form a close working collaboration between the numerical 
relativity and data analysis communities. Several decisions were made that restrict 
the scope of the results reported here: we consider only BBH simulations and have 
not used results from supernova simulations or simulations containing neutron stars; 
the waveform data comes purely from numerical simulations and we do not attempt 
to extend numerical data using post-Newtonian waveforms; the NINJA data set is 
constructed using Gaussian noise to model the response of the Initial LIGO and 
Virgo detectors - no attempt has been made to include non-Gaussian noise transients 
found in real detector data. The comparisons and conclusions reported here are 
thus necessarily limited, and in many cases are only the first steps towards fully 
understanding the sensitivity of data-analysis pipelines to black hole signals. Further 
studies are needed regarding the accuracy and comparison of numerical waveforms, 
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and of how systematic errors in these waveforms can affect parameter estimation. 
Some analyses of numerical waveforms with regard to gravitational-wave detection 
have already been performed [64,72-74], accuracy standards have been developed for 
use of numerical waveforms in data analysis [209] and a detailed comparison of some 
of the waveforms used in the NINJA project was performed in the related Samurai 
project [75]. We expect that subsequent NINJA analyses will build on these results 
to address these issues. 

Despite the limited scope of the first NINJA project, we are able to draw the 
following broad conclusions from this work. Our first conclusion is that the current 
data analysis pipelines used to search LIGO, Virgo and GEO600 data for black hole 
coalescence are able to detect numerical waveforms injected into the NINJA data set 
at the expected sensitivities. Indeed, several of these pipelines are able to detect 
signals that lie outside the parameter space that they target. This is a non-trivial 
statement since most detectability estimates to date for these sources have relied on 
post-Newtonian waveforms, which are valid only when the black holes are sufficiently 
far apart. For many of these pipelines, this is the first time they have been tested 
against numerical waveforms. It should be noted, however, that the NINJA data set 
does not contain non-stationary noise transients so more work is needed to understand 
how detection performance is affected by the noise artifacts seen in real gravitational- 
wave detector data. Our second conclusion is that significant work is required to 
understand and improve the measurement of signal parameters. For instance, among 
the pipelines used in this first NINJA analysis only the Markov-chain Monte-Carlo 
algorithm attempted to estimate the spins of the individual black holes, and the 
estimation of the component masses by the detection pipelines is poor in most cases. 
Improvement in this area will be crucial for bridging the gap between gravitational 
wave observations and astrophysics. NINJA has proven to be extremely valuable at 
framing the questions that need to be answered. 

This paper is organised as follows: In the next section we describe the contributed 
numerical waveforms and Section 3 describes the construction of the simulated 
gravitational-wave detector data used in the NINJA analyses. Descriptions of the 
search methods and results are given in Section 4. The results are grouped by 
search method into search pipelines using modelled waveforms (Section 4.1), search 
pipelines using un-modelled waveforms (Section 4.2) a comparison of inspiral-burst- 
ringdown results (Section 4.3), and Bayesian pipelines (Section 4.4). We conclude 
with a discussion of our results and future directions for NINJA in Section 5. 

2. Numerical Waveforms 

The NINJA project has studied BBH coalescence waveforms submitted by ten 
individuals and teams. Participation in NINJA was open to anyone and the only 
restrictions were that each contribution: (i) was a numerical solution of the full 
Einstein equations, (ii) consisted of only two waveforms, or up to five waveforms 
if they were part of a one-parameter family. No restrictions were placed on the 
accuracy of each waveform. All contributions followed the format specified in [76] . The 
waveforms are plotted in Figures 1 and 2. The contributed waveforms cover a variety of 
physical and numerical parameters. Most simulations model low-eccentricity inspiral, 
the mass ratio q = m\/m,2 ranges from 1 to 4, and the simulations cover a range of 
spin configurations. The initial angular frequency of the I = m = 2 mode ranges from 
0.033/M to 0.203/M (where M denotes the sum of the initial black-hole masses). This 
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initial angular frequency marks where contributors consider the waveform sufficiently 
clean to represent the physical system (e.g. this will be chosen after initial unphysical 
radiation content, often referred to as "junk radiation" in numerical relativity, is 
radiated away). The length of the waveforms varies between a few 100M to over 
4000M. The contributions naturally differ in accuracy, both regarding how well they 
capture the black-hole dynamics and in the extraction of the gravitational- wave signal. 
Table 1 lists a few key parameters that distinguish the waveforms, and introduces 
the following tags for the different contributions and codes: BAM HHB [77-81] 
and BAM FAU [60, 77, 78, 81] are contributions using the BAM code, CCATIE is the 
AEI/LSU code [34,43,52,82,83], Hahndol is the Goddard Space Flight Center's 
code [84,85], LazEv is the RIT code [25,46,86], Lean is Ulrich Sperhake's code [56, 
87,88], MayaKranc is the Georgia Tech/Penn State code [57,74], PU stands for the 
Princeton University code [24,62,89,90], SpEC for the Corncll/Caltcch collaboration 
code [66,91-93], and UIUC stands for the University of Illinois at Urbana-Champaign 
team [94]. 

The codes listed above use different formulations of the Einstein equations, gauge 
conditions, mesh structures, initial data and wave extraction methods; we will attempt 
to give a unified presentation of common features first, and then list further details of 
the approaches separately for each contribution. Full details of each code are given in 
the references. 

The numerical codes follow either of two approaches to solving the Einstein 
equations: (1) the generalised harmonic formulation, which was the basis of Pretorius' 
initial breakthrough simulation of coalescing black holes [24], or (2) the moving- 
puncture approach, following [25,26]. Both approaches result in canonical choices 
for the construction of initial data, the evolution system for the Einstein equations, 
and the treatment of the singularity inside the black- hole horizons. 

2.1. Summary of the simulation algorithms 

2.1.1. Initial Data Due to the presence of constraint equations, specifying initial 
data in numerical relativity is far from trivial, for a general overview see e.g. [101]. All 
of the results presented here make the simplifying assumption of conformal flatness 
for the spatial metric of the initial slice, which leads to some spurious gravitational 
radiation in the initial data. All contributions attempt to model non-eccentric inspiral, 
except for the two data sets PU-T52W and MayaKranc-e02. However, the degree of 
"quasi-circularity" varies, and in general one should bear in mind that the definition of 
eccentricity for fully general-relativistic orbits is not unique (see for example [56, 57]). 
The data set PU-T52W is notable for the fact that the BBH was constructed via 
scalar field collapse. Specifically, the initial data consists of two, compact, dense 
distributions of scalar field energy, separated by some distance and Lorentz boosted in 
opposite directions orthogonal to the line between them. Upon subsequent evolution, 
each scalar field pulse quickly collapses to form a black hole, with all remnant scalar 
field energy radiating away from the domain on the order of the light-crossing time 
of the orbit. This is the same time scale on which spurious gravitational radiation 
present in all current initial-data sets leaves the domain of the inspiral, and hence for 
practical purposes this can be considered a vacuum merger. All other runs start from 
vacuum initial data. 

Most codes (BAM, CCATIE, Hahndol, LazEv, Lean, MayaKranc and the UIUC code) 
adopt the "moving puncture" approach, following [25,26]. These codes use puncture 
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Table 1. Initial conditions for numerical waveforms. The columns list, in 
order from left to right, the name of the contribution or code, the name of the 
run where appropriate, the mass ratio q = m\/m,2 where mi > ma, the spins of 
the black holes in vector form (if only one spin is given, both spins arc equal), 
an estimate of the initial eccentricity of the orbit (the entry qc denotes cases 
where quasi-circular inspiral, i.e. zero eccentricity is modelled, but a value of the 
eccentricity has not been reported), the initial frequency of the (t,m) = (2,2) 
mode (rounded to three digits), the initial coordinate separation of either the 
black-hole punctures or the excision surfaces, and where appropriate the method 
of eccentricity removal. All binaries start out in the xj/-planc with initial momenta 
tangent to the icy-plane. See text for the identification of each contribution, and 
a description of the notation in the last column. The dimensionless spins of the 
BAM FAU run are (-0.634, -0.223,0.333) and (-0.517, -0.542,0.034). 
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Figure 1. Summary of all submitted numerical waveforms: 
r/M Re(h,22) The x-axis shows time in units of M and the j/-axis shows 
the real part of the (£, m) = (2, 2) component of the dimcnsionlcss wave strain 
rh = rfc_|_ — irh x . The top panels show the complete waveforms: the top-left 
panel includes waveforms that last more than about 700M, and the top-right 
panel includes waveforms shorter than about 700M. The bottom panel shows an 
enlargement of the merger phase for all waveforms. 
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Figure 2. Distribution of power into different spherical harmonics. The 

blue line shows (Sf jm |/i£ m r/MI 2 ) 1 ^ 2 . A dashed red line, if present, shows the 
same sum, but excluding the (I, m) = (2, ±2) modes. The separation between 
the two lines gives the relative importance of non (2, ±2) modes. If no red line is 
present for a certain run, then only the (2, ±2) modes were supplied. The layout 
is as in Fig. 1: The top panels show the complete waveforms, whereas the bottom 
panel shows an enlargement of the merger phase. The x-axis shows time in units 
of M. 
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Table 2. Characteristic duration, mass and frequencies of the waveforms 
summarised in table 1. The columns ATioo and /i,ioo gi ve the duration and 
initial frequency of the waveform when scaled to total mass M = 100M© . M30.tr z 
is the total mass of the waveform when it is scaled so that the initial frequency is 
30Hz (this sets the lowest mass at which each waveform can be injected into the 
NINJA data). 

initial data [102-104] to model black holes, resulting in initial data that contain a 
separate asymptotically flat end within each black hole. Constructing such initial 
data is mathematically well understood [103,105]. The codes CCATIE, LazEv, Lean 
and MayaKranc all use the same pseudo-spectral solver for the Einstein constraint 
equations [106], and BAM uses a variant thereof [78]. UlUC-punc initial data is generated 
via the Lorene [107] multi-domain spectral libraries. The Hahndol code uses the 
second-order-accurate multi-grid solver AMRMG [108], which is however tuned to give 
truncation errors typically much smaller than those produced by the evolution code. 

The generalised harmonic codes use conformal thin sandwich initial data [109]. 
PU-CP and SpEC use quasi-equilibrium excision initial data where the interior of 
the black-hole horizons has been excised from the numerical grid. The presence of 
black holes with desired linear momenta and spins is enforced through the boundary 
conditions on the excision surfaces and the numerical outer boundary during the 
solution of the initial-value equations [92,99,110,111]. This "excision technique" is 
based on the defining property of black holes — the horizons act as causal membranes 
and information cannot escape from the inside. The UIUC-cp simulation uses the same 
excised initial data, but fills the BH interior with "smooth junk" , as described in [94], 
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before evolving with the moving puncture technique. 

All codes take input parameters that ultimately determine the individual black- 
hole masses mi, spins Si, momenta Pj and coordinate separation D of the black holes 
(one should however be aware that in the strong field regime of general relativity 
various subtleties are associated with the definition of all of these quantities). In 
addition, the black-hole masses and dimcnsionless spins slowly change during the 
inspiral, which requires additional caution regarding the definition and accuracy of 
the values of mass, spin, etc. There are two common methods to estimate the 
instantaneous individual black- hole masses. One is to calculate the apparent-horizon 
mass, computed from the irreducible mass (given by the area of each hole's horizon) 
and the spin according to Christodoulou's [112] relation mf = m| irr + Sf /(4m? irr ). 
The other, applicable only to puncture data, is to compute the Arnowitt-Deser-Misner 
(ADM) mass [113] at each puncture, which corresponds to spatial infinity in a space 
that contains only that black hole [104]. We generally use the total black-hole mass 
M = TOi+m2 to scale dimensionful quantities, although sometimes the total conserved 
energy (Madm) is used for this purpose. Without loss of generality all codes chose 
the rest frame where Pi = — Pi and, thus, the net linear momentum vanishes initially. 

Those simulations that attempt to model non-eccentric inspiral use initial 
parameters calculated by a number of different methods. Ideal initial parameters 
would produce tangential motion consistent with circular orbits, and radial motion 
consistent with slow inspiral. The various methods to choose initial parameters can 
be broadly characterised as those that attempt to provide only tangential motion (so 
that initially the black holes have no radial momenta), denoted by "T" in the last 
column of Tab. 1, and those that provide both tangential and radial motion (denoted 
by "TR"). The procedures to estimate these parameters are based on properties of 
the initial-data set ( "ID" ) , post-Newtonian methods ( "PN" ) , or an iterative procedure 
following the results of several trial simulations ("it"). In Tab. 1 we indicate which of 
these variants was used, and provide a reference to the specific algorithm; for the post- 
Newtonian methods in particular there are several variants. Note that the estimates 
of the resulting eccentricity range from e ~ 5 x 10 -5 (for the SpEC contribution) up 
to e - 0.02. 

The two data sets from the UIUC contribution actually compare two alternative 
sets of non-spinning, equal-mass, quasi-circular initial data, with initial orbital 
frequency MO = 0.0824: (i) Puncture initial data with coordinate separation 
D/M = 4.369 and initial linear momentum of each BH set according to [100], and (ii) 
Cook-Pfeiffcr initial data with coordinate separation D/M = 4.790 [99, 114] (measured 
from the centroids of the apparent horizons), filling the BH interior with data that 
smoothly connect to the exterior as described in [94]. Both data sets yield the same 
final spin \Sbh \/Mbh = 0.68, but differ at the level of a few percent in radiated 
energy and angular momentum. 

For the eccentric MayaKranc simulation (data set e02), the conservative, third- 
post-Newtonian-order (3PN) expressions in Ref. [115] have been used to specify initial 
data. These expressions require the specification of the eccentricity e and the mean 
motion n = 1ix/T r , where T r is the radial (pericenter to pericenter) orbital period. 
There are three PN eccentricities, which are the same to 1PN order, and we choose 
e t , which appears in the PN Kepler equation, following Ref. [115]. The quantity n 
has been chosen as n = 0.01625/M (T r ~ 387M) and e = 0.2. The binary separation, 
D/M = 15.264, was determined from equation (23) in Ref. [115], and the tangential 
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fa mi n 

0.001M 


BAM HHB 


BSSN 


FD-6 


000 


2 


773 


90 


56, 19 


BAM FAU 


BSSN 


FD-6 


000 


2 


436 


50 


16 


CCATIE 


BSSN 


FD-4 


000 


1 


819 


160 


20 


Hahndol 


BSSN 


FD-4, 6 


000 


2 


> 1000 


45 


19, 13 


LazEv 


BSSN 


FD-4 


ttt 


6 


1281 


40 


3.1 


Lean 


BSSN 


FD-4, 6 


000 


1.25,1 


153.6, 256 


60, 61 


19, 13 


MayaKranc 


BSSN 


FD-4 


000 


2 


317.4 


70 


16, 19 


PU 


GH 


FD-2 


n/a 


n/a 


oo 


50 




SpEC 


GH 


Spectral 


n/a 


n/a 


450 -> 230 


75 - 225 


- 3 


UIUC 


BSSN 


FD-4 


000 


0.25 


409.6 


70 


25 



Table 3. Some properties of the NR evolution codes. The columns list, 
for each contribution, the employed evolution system, the numerical technique 
(FD-k stands for finite differences using k-th order stencils in the bulk), the time 
derivative and r\ choices for the T-driver shift, the approximate location of the 
outer boundary, the radii used for wave extraction, and the finest grid— spacing. 
If two numbers are given they correspond to the two runs of the respective code 
listed in table I (for BAM_HBB, 

^min — 0.019-/V/ applies to all runs with spin). 
For the SpEC run, r max decreases during the run and the waveform is extrapolated 
to r ex t = co based on extraction at radii in the given interval [66, 93]. 

linear momentum, P/M=0.0498, of each black hole at apocenter was obtained from 
J = PD, where J is the total angular momentum computed as a post-Newtonian 
expansion in n and e (equation (21) in Ref. [115]). 

2.1.2. Evolution systems There is a long history of casting the Einstein equations 
into systems of partial differential equations, and in particular into the form of a well- 
posed initial value problem. The process of writing the covariant Einstein equations 
in the form of three-dimensional tensor quantities that evolve in time is commonly 
referred to as a 3+1 split. The fundamental idea is to choose coordinates {x l ,t} 
(i = 1, 2, 3) such that the spacetime metric can be written in the form 

ds 2 = -{a 2 - HjPp^dt 2 + 2j ij f3 j dt dx l + ll0 dx l dx j , (1) 

where 7^ is a positive-definite metric on the slices of constant time t, and the scalar 
function a and vector field (3 l are commonly used to encode the freedom of coordinate 
choice. They may in principle be freely specified, but in practice they are judiciously 
prescribed, usually through further evolution equations. 

The waveforms contributed to NINJA use versions of cither of the two 
formulations for which successful multi-orbit evolutions of black-hole binaries have 
been published so far: the generalised harmonic and the BSSN/moving-puncture 
formulation of the Einstein equations. For overviews of writing the covariant Einstein 
equations as a time evolution problem, see e.g. [116-118]. 

The generalised harmonic formulation (see e.g. [118]) writes the evolution 
equations in manifestly hyperbolic form as a set of coupled wave equations for the 
space-time metric c/ M „. The SpEC code uses this formulation in first order form [119], 
while the PU contribution is based on a second order version of the equations. Gauge 
conditions are enforced by specification of gauge-source functions H^, either as a 
specified function of time, or through evolution equations [24, 66, 93, 120]. 
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All other codes use the first-order-in-time, second-order-in-spacc BSSN 
formulation of the Einstein evolution equations [121-123] in combination with 
hyperbolic evolution equations for the lapse and shift. The BSSN formulation consists 
of making a conformal decomposition of the spatial metric, 7^ = ip 4 Jij , and all other 
variables, and the introduction of T t = dffi , which is treated as an independent 
variable. The moving-puncture treatment of the BSSN system involves evolving 
not the conformal factor ip but cither cj) = hiip (CCATIE), W = ip~ 2 (BAM FAU, 
Hahndol [55, 124]), or \ — 4>~ A (used by all other BSSN codes); it also consists of the 
gauge choices that we will summarise next. 

All BSSN-based contributions evolve the lapse according to the 1+log slicing 
condition [125], 

(St -ftdi)a. = -2aK. (2) 

The shift vector field (3 l is evolved according to some variant of the T-driver condition 
[83,85]). During the evolution these gauge conditions change the geometry of the 
"puncture singularity" and soften the singularity as discussed in [126-129]. 
The original T-driver condition introduced in [83] is 

d t p = -B\ d t B l =d t r - V B\ (3) 

The factor of 3/4 is chosen such that at large distances the propagation speed of the 
hyperbolic equation (3) equals the coordinate speed of light [83] , and the quantity n is 
a parameter with the dimensions of the inverse of a mass and affects coordinate drifts: 
larger values of 77 lead to a stronger initial growth of the apparent horizon, and thus 
to a magnification effect for the black holes [77]. Variants of this condition [25,26, 
85, 130, 131] consist of replacing some or all of the d t derivatives with do = d t ~ (i l di. 
We will label these options with reference to each of the three time derivatives in (3): 
"ttt" denotes that dt is used for all three derivatives, "000" denotes usage of do- The 
properties of the different choices are studied in [85, 131], and in [131] it is proven that 
the combination of the BSSN equations with the "1+log" slicing condition (2) and 
the "000" shift choice yields a well-posed initial-value problem. 

Small differences in the evolutions also originate in the choice of initial lapse 
(all BSSN codes initialise the shift quantities (3 l and B % to zero). We first define a 
Brill-Lindquist-likc conformal factor, ipBL = 1 + rn\ >v /2r\ +m2 jP /2r2, where is the 
distance to the Ath. puncture, and m\ tP and m2. p parametrise the masses of the black 
holes, although they are not in general equal to mi and m-i- The RIT contributions 
choose a(t = 0) = 2/(1 + ip% L ), as does the Hahndol-non contribution, while the 
Hahndol-kick contribution uses an approximate a(t = 0) derived from the late-time 
"1+log" Schwarzschild slicing [127]. BAM HHB, MayaKranc and the UIUC group use 
a(t = 0) = ipg 2 L , and BAM FAU choose a(t = 0) = [{ip B L - l)/2 + l]" 4 . 

The generalised harmonic codes (PU and SpEC) employ black- hole excision, i.e., 
they excise from the computational grid a region around the singularities inside each 
black hole. 

2.1.3. Radiation Extraction All groups use one of two popular methods to estimate 
the gravitational- wave signal at a finite distance from the source: The SpEC and CCATIE 
contributions use the Zcrilli-Moncricf/Sarbach-Tiglio pcrturbative formalism [132- 
134] (with SpEC following a version restricted to a Minkowski background in standard 
coordinates [135]), all other contributions use the Newman- Penrose curvature scalar 
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Tp4. Both methods are implemented in the CCATIE code, and have been shown to 
give similar results [34,43]). Summaries and details on the implementations within 
particular codes can be found, for instance in the references listed in Table 1. Since 
the gravitational-wave signal can only be defined unambiguously at null infinity, one 
typically considers several extraction radii and performs some form of convergence 
test, although for the present purpose most groups only report results for a single 
extraction radius. At finite radius both methods depend on the coordinate gauge, 
and the Newman-Penrose method additionally requires the choice of a tetrad, which 
is obtained by Gram-Schmidt orthonormalisation of a tetrad of coordinate vectors. 

For this work, all waveforms have been contributed as spherical harmonic modes 
of spin- weight —2 of the strain, according to the specification in [76]. Computation 
of the strain from the Zerilli-Moncrief odd- and even-parity multipoles of the metric 
perturbation requires one time integration [43,133], in the Sarbach-Tiglio formalism 
the strain is algebraically related to the invariants at leading order in the inverse 
radius [133, 136], and computation of the strain from the Newman-Penrose curvature 
scalar i/>4 requires two time integrations. Time integration requires the proper 
choice of integration constants, and may require further "cleaning procedures" to 
get rid of artifacts resulting from the finite extraction radii. For example, for the 
BAM HHB contribution unphysical linear drifts were removed by a variant of the 
method described in [70], where higher order than linear polynomials were used to 
remove unphysical drifts from higher modes to further improve the properties of the 
derived strain. In the RIT contribution, the strain was computed by taking the Fourier 
transform of ip±, removing modes in a small region around u> = 0, then dividing by 
— uj 2 and taking the inverse Fourier transform. 

2.1. 4-- Numerical Methods and Computational Infrastructure There arc large overlaps 
regarding the numerical methods in the present waveform contributions. With the 
exception of the SpEC code, which uses a multi-domain pseudo-spectral method, all 
codes use finite-difference methods to discretise the equations. With the exception of 
the PU contribution, which uses a second-ordcr-accuratc implicit evolution scheme, 
all other codes use an explicit algorithm based on method of lines: Usually standard 
fourth-order-accurate Runge-Kutta time stepping, except for the SpEC code which uses 
a fifth order Cash-Karp time-stepper with adaptive step-size. 

The moving-puncture/BSSN-based codes use standard centred finite differencing 
stencils; however the terms corresponding to the Lie-derivative with respect to the 
shift vector are off-centred (up-winded) by one grid-point. The CCATIE, MayaKranc, 
LazEv and UIUC codes use fourth-order-accurate stencils, the BAM code uses sixth- 
order stencils, the Halmdol code uses sixth-order stencils combined with fifth-order 
up- winded stencils [137], and the Lean code uses fourth-order for equal-mass and 
sixth-order for unequal-mass data sets. All of these codes add standard fifth-order 
Krciss-Oliger dissipation [138, 139] to the right-hand-sides of the evolution equations. 
The finite-difference orders described here apply to the bulk of the computational 
domain. There are contributions at other orders in different parts of the codes, which 
we will describe below. However, the finite-difference order in the bulk plays the 
dominant role in defining the accuracy of the present simulations (and indeed the 
spatial finite-differencing order seems to dominate over the order of time integration 
when sufficiently small time steps are used), and for that reason we list in Tab. 3 the 
bulk spatial finite-difference order. 

All codes except the SpEC code use variants of Bcrgcr-Oligcr mesh- refinement. 
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The PU and Hahndol codes employ full adaptive mesh refinement, while the other 
codes use a hierarchy of fixed refinement boxes which follow the motion of the black 
holes. Several of the codes are based on the Cactus computational toolkit [140, 
141] and the Carpet mesh-refinement code [142,143] (CCATIE, Lean, MayaKranc, 
LazEv, UIUC). The BAM HHB and BAM FAU contributions both use the BAM mesh 
refinement code. The Hahndol code uses the PARAMESH infrastructure [144] with a 
uniform time step; all other mesh refinement codes use a time step that depends on 
the grid spacing, and for these codes time interpolation at mesh-refinement boundaries 
introduces second-order errors. 

For interpolation between meshes of different spacing, the groups that used fourth- 
or higher-order methods all use fifth-order-accurate (CCATIE, UIUC, LazEv, Lean, 
MayaKranc and Hahndol's 4:1 "non" data) or sixth-order-accurate (BAM and Hahndol's 
3:1 "kick" data) polynomial interpolation in space between different refinement levels 
so that all spatial operations of the AMR method (i.e., restriction and prolongation) 
are sixth-order accurate and the second derivatives of interpolated values are at least 
fourth-order accurate. 

A proper numerical treatment of gravitational waves in asymptotically flat 
spacetimes would include null infinity and not require boundary conditions at some 
finite distance from the source. Most codes circumvent this problem in essentially 
heuristic ways. The PU code uses spatial compactification combined with numerical 
dissipation, all BSSN codes use heuristic outgoing wave boundary conditions (which 
will in general violate constraint preservation and potentially well-posedness and 
will result in reflections of the outgoing radiation). The SpEC code, in contrast, 
uses constraint-preserving outer boundary conditions which are nearly transparent 
to outgoing gravitational radiation and gauge modes [145]. 

Note that several of the groups use the same apparent horizon finder code 
(AHFlNDERDlRECT) [146,147] (Hahndol, UIUC, CCATIE, LazEv, MayaKranc, Lean). 

2.2. Accuracy 

Estimates on accuracy are reported for the BAM HHB and SpEC contributions. For 
the BAM HHB simulations reasonably clean sixth-order convergence was observed, as 
reported in [79,80]. In the waveform rVl^, extracted at R ex = 90M, the uncertainty 
due to numerical errors and the use of finite extraction radii is estimated as 0.25 
radians in the phase and less than 3% in the amplitude of the / = 2, m = 2 mode. 
Modes up to I — 8 were calculated; the relative phase uncertainty is the same for all 
of them (the absolute phase uncertainty is proportional to m), but we estimate that 
the amplitude uncertainty increases to as much as 10% for the highest modes. The 
SpEC contribution is the only one that extrapolates the gravitational wave signal to 
infinite extraction radius (using third-order polynomial extrapolation [66]). Various 
convergence tests indicate that the resulting extrapolated waveform is accurate to 0.02 
radians in phase and 0.5 percent in amplitude [66]. 

3. Construction of the NINJA data set 

The data provided by the numerical relativity groups follows the format outlined 
in [76] , which is based on the mode decomposition of the gravitational radiation field 
at large distances from the source. If we specify a gravitational waveform in the 
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Transverse- Traceless (TT) gauge, we only need the spatial components fty . We assume 
that we are sufficiently far away from the source so that the 1/r piece dominates: 

^•=^ + 0(r- 2 ) , (4) 

where M is the total mass of the system, r is the distance from the source, and Aij is 
a time-dependent TT tensor. In the TT gauge, hij has two independent polarisations 
denoted h + and h x and the complex function h + — ih x can be decomposed into modes 
using spin-weighted spherical harmonics ~ 2 Yi m of weight -2: 

, » oo £ 

h+-ih x = —J2 H ^ ~ 2Y *™(h <f>) ■ (5) 

e=2 m=-l 

The expansion parameters Hi m are complex functions of the retarded time t — r, 
however if we fix r to be the radius of the sphere at which we extract waves then Hi m 
are functions of t only. The angles i and <j> are respectively the polar and azimuthal 
angles in a suitable coordinate system centred on the source. This decomposition 
is directly applicable to non-precessing binaries. Otherwise, a comparison of the 
waveforms requires a careful treatment of mode-mixing effects due to rotations of 
the frame; see for instance [148]. The numerical data contributed to NINJA is given 
in the form of an ASCII data file for each mode (t, m), with accompanying meta-data 
describing the simulation [76]. Only modes that contribute appreciably to the final 
waveform are included, at the discretion of the contributing group. Each data file 
consists of three columns: time in units of the total mass, and the real and imaginary 
parts of the mode coefficients H( m as a function of time. Note that the total mass M 
scales both the time and the amplitude; thus the BBH waveforms for each simulation 
can be scaled to an arbitrary value of the mass. (This is not true in the case of 
simulations which include matter fields, but we do not consider such waveforms here.) 

To model the signal seen by a gravitational-wave detector, we need to calculate 
the detector strain h(t) from the above mode decomposition. To do this, we must 
choose particular values of the total mass, orientation and distance from the detector. 
Given the H(, m , the total mass, the distance to the source, and the angles (i, <j>), we 
calculate h + ^ x using Equation. (5), and use the detector response functions i*+ )X (see, 
for example, Ref. [1]) to calculate the observed strain 

h(t) = h+{t)F+{a, 6, VO + h x (t)F x (q, S, t/>) . (6) 

Here (a, 6) are sky-angles in the detector frame, ip is the polarisation angle and the 
time t is measured in seconds. In this analysis, we wish to simulate signals that might 
be observed by the Initial LIGO and Virgo detectors. There are three LIGO detectors: 
a 4 km detector and a 2 km detector at the LIGO Hanford Observatory (called HI and 
H2, respectively) and a 4 km detector at the LIGO Livingston Observatory (called LI). 
The Virgo detector is a 3 km detector in Cascina, Italy (called VI). We used the same 
two-letter codes for the simulated NINJA detectors. Since the location and alignment 
of the three observatories differ, we must use the appropriate detector response and 
arrival time to compute the strain waveform h(t) seen at each observatory. This 
ensures that the waveforms arc coherent between the detectors and simulate a true 
signal. 

To model the detector noise, we generated independent Gaussian noise time scries 
n(t), sampled at 4096 Hz, for each detector. This sample rate was chosen to mimic 
that used in LSC- Virgo searches and assures a tolerable loss in signal-to-noise ratio 
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Frequency [Hz] 

Figure 3. The NINJA data noise curves and the design spectra of the first 
generation LIGO and Virgo detectors. 

due to the discrete time steps. Stationary white noise time series are generated and 
coloured by a number of time-domain filters designed to mimic the design response 
of each of the LIGO and Virgo detectors. Figure 3 shows the one-sided amplitude 
spectral density yj S n (f) of each time detector's time series, where S n (f) is defined by 

(nimn = \sn{\fMf - n- (?) 

n(/) denotes the Fourier transform of n(t) and angle brackets denote averaging over 
many realisations of the noise. We see from Figure 3 that the noise power spectrum 
of the NINJA data set closely approximates the Initial LIGO design sensitivity in the 
frequency range of interest (30-10 3 Hz). There is a slight discrepancy with the Virgo 
design curve at low frequencies (between approximately 20 and 150 Hz), which is an 
artefact of the Virgo noise generation procedure. Narrow-band features such as the 
violin and mirror modes were removed from the detector response used to compute 
the NINJA data, but were included in the calculation of the Virgo design curve [208] . 
The 1/ / tails of these narrow-band features are responsible for the small discrepancy. 

Having generated the simulated detector data, we then generated a population of 
simulated signals using the numerical relativity data. This population was constructed 
to cover a broad range of masses and signal amplitudes. We required that the starting 
frequency of the dominant I = m = 2 mode of the signal was not more than 30 Hz, 
an appropriate threshold given the sensitivity curve of the Initial LIGO and Virgo 
detectors. This sets a minimum mass at which each waveform can be injected, 
which is given in Table 2. The minimum possible injection mass is therefore 36M Q . 
The maximum mass was chosen as 350 M Q . To get a good sample of long injected 
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waveforms, we systematically chose a lower range of masses for the longer waveforms. 
No restrictions were placed on the other simulation parameters, i.e., the spins, mass- 
ratios and eccentricities. We ensured that waveforms from all the participating groups 
were equitably represented by generating approximately 12 signals from the waveforms 
supplied by each group. The time interval between adjacent injected signals was chosen 
to be a random number in the range 700 ± 100 s. 

Given these constraints, we generated the parameters of the signal population. 
The logarithm of the distance to the binary was drawn from a uniform distribution 
ranging from 50 Mpc to 500 Mpc, and the source locations and orientations were drawn 
from an isotropic distribution of angles. We then computed waveforms corresponding 
to this population and at the appropriate sampling rate. We required that the optimal 
matched filter signal-to-noise ratio of any injection be greater than five in at least one 
of the four simulated detectors. Any waveform that did not satisfy this constraint was 
discarded from the population. Subject to this condition, the distances of injected 
signals varied from 52 Mpc to 480 Mpc (median at 145 Mpc), the injected total mass 
range was 36M Q < M < 346M Q (median at 155M©), with individual component 
masses in the range 11M Q < nij < 193A-/ Q . 

Finally, the waveforms h{t) were added to the simulated detector noise n(t) to 
generate the NINJA data set s(t) — n(t) + h(t). As described above, care was taken 
to ensure that signals were coherently injected in the data streams from the four 
detectors. The software for carrying out this procedure is freely available as part of 
the LSC Algorithm Library (LAL) [149]. 

The data set used in this analysis consisted of a total of 126 signals injected in a 
total of 106 contiguous segments of noise each 1024 s long, thus spanning a duration 
of a little over 30 hours. Figure 4 shows the mass, spin and distance of the waveforms 
contained in the NINJA data set. 

4. Data Analysis Results 

Analysis of the NINJA data was open to all and nine groups submitted contributions 
using a variety of analysis techniques. Participating groups were provided with 
the NINJA data set containing signals embedded in noise and the parameters of 
the injected signals. Analysts were not given access to the raw numerical-relativity 
waveforms or noiseless injection data. 

Methods used to analyse the NINJA data include: matched-filter based searches, 
un-modelled waveform searches using excess-power techniques, and Bayesian model- 
selection and parameter-estimation techniques. Where possible, the performance of 
different searches is compared. The limited scope of the NINJA data set makes detailed 
comparisons difficult, however. A list of the data-analysis contributions is shown in 
Table 4. 

In sections 4.1 and 4.2 we describe results of analyses using modelled (matched- 
filter) and un- modelled waveforms, respectively. Comparisons between these analyses 
are given in Section 4.3 and Section 4.4 presents the results of Bayesian model-selection 
and parameter-estimation analyses. 

4-1- Search pipelines using modelled waveforms 

When the waveform of the target signal is known, matched filtering is the optimal 
search technique for recovering signals buried in stationary noise [150,151]. This 
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Figure 4. The total mass and distance of the 126 NINJA injections. 

The grey scale encodes the sum of the dimensionless spins of the black holes, 
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Group 


Analysis 


Section 


AEI 


Phcnomcnological Waveforms in CBC pipeline 


4.1.3 


Birmingham 


Baycsian Model Selection 


4.4.2 


Cardiff 


Post-Newtonian (PN) Templates in CBC pipeline 


4.1.1 


Cardiff. Maryland 


EOBNR waveforms in CBC pipeline 


4.1.3 


Goddard 


Hilbcrt Huang Transform 


4.2.2 


Northwestern 


Markov Chain Monte Carlo 


4.4.1 


Syracuse 


Extended rj PN Templates in CBC pipeline 


4.1.1 


UMass, Urbino 


Q-pipeline analysis 


4.2.1 


UWM 


PN templates in CBC pipeline, Neyman-Pearson criteria 


4.1.2 


UWM, UMass, Urbino 


Ringdown analysis 


4.1.4 


UWM, UMass, Urbino 


Inspiral, Merger, Ringdown combined search 


4.3 



Table 4. The data-analysis contributions to the NINJA project. "CBC 
pipeline" refers to the LSC- Virgo Compact Binary Coalescence group's analysis 
pipeline, described in section 4.1. 
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section describes the results of filtering the NINJA data with matched-filter based 
analysis pipelines. Results are given for waveforms that span only the inspiral 
signal, the ringdown alone, and the full inspiral, merger and ringdown. Although the 
morphologies of these waveforms differ, the underlying analysis techniques are similar 
in all cases. All the contributions in this section use a pipeline developed by the LSC 
and Virgo Collaboration to search for gravitational waves from binary neutron stars 
and black holes in a network of detectors [17,152]. We first describe the features of 
this pipeline common to all the contributed matched-filter analyses before presenting 
the results of searching the NINJA data using different matched-filter templates. 

The LSC- Virgo search pipeline performs a series of hierarchical operations in 
order to search for real signals buried in the detector noise: Given a desired search 
parameter space and waveform model, a "bank" of templates is created to cover the 
parameter space such that the fractional loss in signal-to-noise ratio (SNR) between 
any signal and the nearest template is less than a specified value (typically 3%). All 
the NINJA inspiral searches use a non-spinning template bank parametrised by the 
two component masses of the binary [153-155]. It has been found that inspiral searches 
for spinning binaries using waveforms which neglect the effect of spin are reasonably 
effective in most cases [152,156]. Ringdown searches use a two parameter template 
bank parametrised by the frequency and quality factor of the signal constructed to 
cover the desired range of mass and spin [157]. Data from each of the detectors is 
separately match filtered against this bank of waveforms [157, 158] and a "trigger" is 
produced whenever the SNR exceeds the desired threshold. All the analyses used a 
threshold of 5.5. A test is then performed which discards triggers which do not have 
coincident parameters in two or more detectors (time and masses for inspiral searches, 
and time, mass and spin for ringdown searches) [159, 160]. These coincident triggers 
provide the gravitational- wave candidates for the ringdown analysis. The triggers are 
ranked by a detection statistic p c constructed from the SNRs of the N > 2 individual 
triggers in a coincidence by p c = (YliLi Pi) 1 ^ 2 - Coincident inspiral triggers are subject 
to a second stage of filtering in which "signal-based vetoes" are also calculated, which 
aim to separate true signals from noise fluctuations. These include the \ 2 [161] and 
r 2 [162] tests. Signal-based vetoes could also be employed for ringdown searches, but 
at present they are not implemented in the pipeline. For each trigger, we construct 
an effective SNR p c g, which combines the matched-filter SNR and the value of the \ 2 
signal based veto [161]. Explicitly, the effective SNR is defined as [17, 152] 



where DOF signifies the number of degrees of freedom in the x 2 test. For signals 
of moderate SNR, which are a good match to the template waveform, the expected 
value of the \ 2 is unity per degree of freedom and consequently the effective SNR 
is approximately equal to the SNR. Non-stationarities in the data typically have 
large values of \ 2 an d consequently the effective SNR is significantly lower than 
the SNR. A second test is then performed to discard coincidences in which signal- 
based vetoes reduce the number of triggers to less than two. These coincidences 
provide the candidate gravitational wave signals for the inspiral-bascd pipelines and 
they are ranked by the combined effective SNR p c g = (X^ili /°eff J 1 ^ 2 - To evaluate 
the sensitivity of the analyses, we compare the list of gravitational-wave candidates 
generated by filtering the NINJA data to the parameters of the inject numerical 
relativity signals. 




(8) 
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Six groups contributed matched-filter results to this analysis and the results can 
be roughly divided into three categories based on the waveform templates used: (i) 
searches based on the stationary-phase approximation to the inspiral signal, which are 
designed to capture various stages of the inspiral, merger and ringdown, (ii) searches 
which use waveforms designed to model the full inspiral-mcrgcr-ringdown signal, (iii) 
searches using ringdown-only waveforms obtained from black hole perturbation theory. 
Within these categories, different parameter choices were made in order to investigate 
the ability of the pipeline to detect the numerical relativity simulations. Each of these 
three approaches is described independently in the following sections. A comparison 
between these results is given in Section 4.3. 

If.. 1.1. Stationary Phase Inspiral Templates The workhorse template of the LSC- 
Virgo search pipeline is based on the stationary-phase approximation to the Fourier 
transform of the non-spinning post-Newtonian inspiral [158,163]. This waveform 
(referred to as SPA or TaylorF2) has been used in the search for binary neutron 
stars [13-15, 17] , sub-solar mass black holes [13, 16, 17] and stellar mass black holes [13] . 
The TaylorF2 waveform is parametrised by the binary's component masses m\ and 
m 2 (or equivalently the total mass M = mi + ?ti 2 and the symmetric mass ratio 
77 = m 1 m 2 /M 2 ) and an upper frequency cutoff f c . Amplitude evolution is modelled 
to leading order and phase evolution is modelled to a specified post-Newtonian order. 
In this section we investigate the performance of TaylorF2-based searches on the three 
simulated LIGO detectors. Results which include the simulated Virgo detector are 
described in the next section. Several analyses were performed which test the ability 
of TaylorF2 waveforms to detect numerical relativity signals. The analyses differed in 
the way the TaylorF2 waveforms or the template bank were constructed. The results 
of these searches are summarised in Table 5, each column giving the results from 
a different search with a summary of the chosen parameters. We first describe the 
parameters varied between these analyses and then present a more detailed discussion 
of the results. 

All TaylorF2 NINJA analyses used restricted templates (i.e. the amplitude is 
calculated to leading order), however the phase was calculated to various different post- 
Newtonian orders [164]. Phases were computed to either two [165, 166] or three point 
five post-Newtonian order [167-169] since these are, respectively, the order currently 
used in LSC- Virgo searches [13] and the highest order at which post-Newtonian 
corrections are known. After choosing a post-Newtonian order, one chooses a region of 
mass-parameter space to cover with the template bank. Figure 5 shows the boundaries 
of the template banks used in the analyses. One search used the range used by 
the LSC- Virgo "low-mass" search [13] (mi, 777,2 > 1M ,M < 35M ) and all other 
searches used templates with total masses in the range 20M Q < M < 90M Q . These 
boundaries were chosen since there were no signals in the NINJA data with mass 
smaller than 36M Q and there is little, if any, inspiral power in the sensitive band of 
the NINJA data for signals with M > 1OOM . The standard LSC- Virgo template 
bank generation code [154] restricts template generation to signals with 77 < 0.25, 
since it is not possible to invert M and r\ to obtain real-valued component masses for 
77 > 0.25. All but one of the searches enforced this constraint, with the 0.03 < 77 < 0.25 
for the low-mass CBC search and 0.1 < 77 < 0.25 for the other "physical-77" searches. 
It is, however, possible to generate TaylorF2 waveforms with "unphysical" values of 
77 > 0.25. In two separate studies using Goddard and Pretorius waveforms [64], and 
Caltech-Corncll waveforms [72] it was observed that match between numerical signals 
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Analysis 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


Freq. Cutoff 


ISCO 


ISCO 


ERD 


ERD 


WRD 


WRD 


PN Order 


2 PN 


2 PN 


2 PN 


3.5 PN 


3.5 PN 


3.5 PN 


Total Mass M 


2-35 


20-90 


20-90 


20-90 


20-90 


20-90 


i] range 


n HQ n or; 
U.Uo-U.zo 


U.lU-U.zo 


n i n nor; 
U.lU-U.zo 


n in n o K 
U.lU-U.zo 


(i in n 

U.lU-U.zo 


n in i 

U.1U— 1 


Found Single 
(HI, H2, LI) 


69, 66, 75 


72, 43, 66 


83, 51, 81 


91, 56, 87 


90, 55, 88 


90, 56, 88 


Found 
Coincidence 


49 


59 


79 


82 


82 


84 


Found Second 
Coincidence 


48 


59 


77 


81 


81 


81 



Table 5. Results of inspiral searches using TaylorF2 templates. There 
were 126 injections performed into the data. The table above shows the number of 
injections which were recovered from the three simulated LIGO detectors (HI, H2 
and LI) using various different waveform families, termination frequencies /iscOi 
/erd and /wrd (as described in the text), and post-Newtonian orders. 



and TaylorF2 templates could be increased by relaxing the condition r\ < 0.25. One 
NINJA contribution uses a template bank with 0.1 < ?y < 1.0 to explore this. 

Finally, it is necessary to specify a cutoff frequency at which to terminate the 
TaylorF2 waveform. In the LSC- Virgo analyses, this is chosen to be the innermost 
stable circular orbit (ISCO) frequency for a test mass in a Schwarzschild spacetime 

c 3 

/lSC ° = dTd^J- (9) 
This cutoff was chosen as the point beyond which the TaylorF2 waveforms diverge 
significantly from the true evolution of the binary [164]. More recently, comparisons 
with numerical relativity waveforms have shown that extending the waveforms up 
to higher frequencies improves the sensitivity of TaylorF2 templates to higher mass 
signals [64,72]. The NINJA TaylorF2 analyses use templates terminated at the 
ISCO frequency and two additional cut-off frequencies: the effective ringdown (ERD) 
frequency and a weighted ringdown ending (WRD) frequency. The ERD frequency 
was obtained by comparing post-Newtonian models to the Prctorius and Goddard 
waveforms [64] . The ERD almost coincides with the fundamental quasi-normal mode 
frequency of the black hole formed by the merger of an equal-mass non-spinning black- 
hole binary. The weighted ringdown ending (WRD) frequency lies between ISCO and 
ERD, and was obtained by comparing TaylorF2 waveforms to the Caltech-Cornell 
numerical signals [72] . 

The results of these searches are reported in Table 5. The principal result is 
the number of injected signals detected by the search. For simplicity, we define 
a detected signal as one for which there is a candidate gravitational-wave signal 
observed within 50 ms of the coalescence time of the injection, determined by the 
maximum gravitational-wave strain of the injected signal. We do not impose any 
additional threshold on the measured SNR or effective SNR of the candidate. For a 
single detector, this will lead to a small number of falsely identified injections, but for 
coincidence results the false alarm rate is so low that we can be confident that the 
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Figure 5. Boundaries of the template banks used in inspiral searches 

as a function of total mass M and symmetric mass ratio 77. The crosses show 
the location of the injections in the NINJA data set. The numbers in the legend 
correspond to entries in table 5. Bank 6 extends in a rectangle up to r] = 1.00, 
as indicated by the arrows. NP is the bank used in the Neyman-Pearson analysis 
described in Section 4.1.2. 



triggers are associated with the injection. We now describe these results in the order 
that they appear in Table 5. 

Search (1) used second order post-Newtonian templates terminated at /isco with 
a maximum mass of M < 35M Q . Despite the fact that no NINJA injections had a 
mass within the range of this search, a significant number of signals were still recovered 
in coincidence both before and after signal consistency tests. Although the templates 
are not a particularly good match to the injected signals, they are still similar enough 
to produce triggers at the time of the injections. Search (2) changed the boundary of 
the template bank to 20Mq < M < 90Mq, but left all other parameters unchanged. 
The number of detected signals increases significantly as more signals now lie within 
the mass range searched. 

Search (3) extended the upper cutoff frequency of the waveforms to /erd- 
The number of signals detected increased from 59 to 77, as expected since these 
waveforms can detect some of the power contained in the late inspiral or early 
merger part of the signal [64,72]. Search (4) extends the post-Newtonian order to 
3.5 PN, slightly increasing the number of detected signals to 81. With the limited 
number of simulations performed in this first NINJA analysis, it is difficult to draw 
a strong conclusion, although there does seem to be evidence that the higher post- 
Newtonian order waveforms perform better, consistent with previous comparisons of 
post-Newtonian and numerical relativity waveforms [64,71,72,79,170]. Search (5) 
uses an upper-frequency cutoff of /wrd for the templates. The number of injections 
found in coincidence for this search is the same as the search using 3.5 order templates 
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Figure 6. Results from the extended template bank. Left: The template 
bank generated by the LSC- Virgo search pipeline (circles) and the bank obtained 
by extending to t) < 1.00 (crosses). In this figure the bank is parametrised 
by to and T3 which are related to the binary masses by r = 5M/(256f?^) 
and T3 = irM/(8r]v^ ), where vo = (ttM/o) 1 ' 3 is a fiducial velocity parameter 
corresponding to a fiducial frequency /o = 40.0Hz. Right: The signal-to-noise 
(SNR) ratio at which NINJA injections were recovered using the 77 < 0.25 bank 
(squares) and the rj < 1 extended bank (circles) in the Hanford detectors, given 
by P = (Phi P^) 1 ' 2, The SNR of the signal recovered using the extended 
bank shows with significant (> 10%) increases over the standard bank for certain 
injections. 



with a cutoff of /erd; although there are slight differences in the number of found 
injections at the single detector level. 

Search (6) extends the template bank of search (5) to unphysical values of the 
symmetric mass ratio. Extending the bank to rj < 1 increases the number of templates 
in the bank by a factor of ~ 2. The original and modified template banks are shown 
in Figure 6. With the extended template bank the number of injections found in 
coincidence remains the same as search (5) after signal-based vetoes are applied. 
However, many of the injections are recovered at a higher SNR, particular the low- 
mass signals, as shown in Figure 6. Some injections show a reduction in SNR; more 
work is needed to understand this effect. 

Finally, we note that the majority of signals passed the \ 2 signal-based veto with 
the thresholds used in the LSC- Virgo pipeline. The last two lines of Table 5 show the 
number of recovered signals before and after these signal-based vetoes are performed. 
The post-Newtonian templates and numerical relativity signals are similar enough 
that virtually all of the injected signals survive the signal based vetoes. 

To illustrate the results of these analyses in more detail, Figure 7 shows which 
signals were detected and which were missed by the 3.5 order post- Newtonian TaylorF2 
templates terminated at /erd, as a function of injected total mass and effective 
distance of the binary (a measure of the amplitude of the signal in the detector), 
defined by [158] 

D cB = d I yjFl{l + cos 2 t) 2 / 4 + Fl cos 2 1 , (10) 

where d is the luminosity distance of the binary. 

One signal, with total mass of 110M Q and effective distance ~ 200 Mpc, was 
missed while others with similar parameters were found. This signal was one of the 
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Figure 7. Found and missed injections using TaylorF2 templates 
terminated at ERD, plotted as a function of the injected effective distance 
in Hanford (left) and Livingston (right) and the total mass of the injection. Since 
the LIGO Observatories are not exactly aligned, the effective distance of a signal 
can differ, depending on the sky location of the signal. The vertical bars mark 
the limits of the template bank used in the search. For the lower masses, we see 
that the majority of the closer injections are found in coincidence in all three of 
the detectors. There is then a band of injections which are found only in two 
detectors - Hf and LI and not the less sensitive H2 detector. For higher masses, 
the results are less meaningful as the template bank was only taken to a total 
mass of 9OM . 



Princeton waveforms (labelled PU-eO . 5 in Figure 1) for which the maximum amplitude 
occurs at the start of the waveform rather than at coalescencej, rendering our simple 
coincidence test invalid. The injection finding algorithm compares the peak time to 
the trigger time and, even though triggers are found at the time of the simulation, 
there are no triggers within the 50 ms window used to locate detected signals. 

Figure 8 shows the accuracy with which the total mass and coalescence time of the 
binary are recovered when using the 3.5 post-Newtonian order Taylor F2 templates. 
The total mass fraction difference is computed as (Mjnjected - -^detected) /-^injected- For 
lower mass signals, the end time is recovered reasonably accurately, with accuracy 
decreasing for the high mass systems. The total mass recovery is poor for the 
majority of signals, with good parameter estimation for only a few of the lowest mass 
simulations. 

4-1-2. Four-detector Inspiral Search The inspiral analysis described in Section 4.1.1 
considered data from the three simulated LIGO detectors. We now extend the 
analysis to include data from the simulated Virgo detector. In addition, we impose an 
alternative criterion, based on the Neyman-Pearson formalism [151], to determine 
those injections which were detected by the pipeline. In the previous section an 
injection was classified as found by the search if gravitational-wave candidate existed 
within 50 ms of the peak time of the numerical data. Here, we consider a signal to be 
found is there is an associated candidate whose significance exceeds a pre-determined 

| That the maximum occurs at the start of the waveform is in part an "artifact" of the double-time 
integration from the Newman-Penrose scalar ip4 to the metric perturbation h, and in part a coordinate 
artifact. The two integration constants were chosen to remove a constant and linear-in-time piece 
for h, however, there is still a non-negligible quadratic component; we suspect this is purely gauge, 
though lacking a better understanding of this it was not removed from the waveform. 
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Figure 8. Parameter accuracy using TaylorF2 templates terminated at 
ERD.Left: Accuracy with which the total mass is recovered. The template 
bank covers the region 20Mq < M < 90Mq, hence the mass of injections 
with M > 90Mq are always underestimated. Even within the region covered 
by the bank, the TaylorF2 templates systematically underestimate the mass of 
the injected signals and the total mass is recovered accurately only for a few 
injections. The vast majority of recoverd signals have an error of 40% or greater. 
Right: Accuracy of determining the coalescence time of the injections. The end 
time is not recovered accurately, the timing error can become as large as 50ms, 
the limits of the injection window. 



threshold. Specifically, we require the candidate to have a significance greater than 
any candidate arising due to noise alone. This allows us to probe in more detail the 
effect of signal-based vetoes and the efficaciousness of the effective SNR statistic in 
analysis of the NINJA data. 

Data from all four simulated NINJA detectors was analysed using the CBC 
pipeline as described in column (1) of Table 5. In addition, a second analysis 
was performed with that the template bank extended to cover the region from 
2Mq < M < IOOMq, with all other parameters unchanged. The search can therefore 
be though of as the simplest extension of the standard LSC- Virgo "low mass" CBC 
search [13]. The boundary of the template bank used is shown in Figure 5. 

In this analysis, we choose a detection statistic and claim that a gravitational-wave 
candidate is present if the value of this statistic exceeds a pre-dctermined threshold. 
All candidates are considered detections. The threshold is chosen so that the false 
alarm probability — the probability that a noise event will be mistaken for a real 
signal — is tolerable. The efficiency of this method depends on how close the chosen 
statistic is to the optimal detection statistic. It is well known that the matched 
filter SNR is the optimal statistic for known signals in a single detector if the noise 
is stationary [150,151]. For a network of detectors containing stationary noise, the 
optimal statistic is the coherent signal-to-noise ratio Pcoherent [171]. At the time of 
this analysis, calculation of /^coherent was not available in the CBC pipeline, so we 
instead compute a combined SNR from the i detectors, p c — (X^ili PiY^ 2 > as a simple 
alternative. In the presence of non-Gaussian noise, the effective SNR, described in 
Section 4.1, has shown to be an effective detection statistic [17]. In this analysis, we 
also consider the combined effective SNR p c s = Q3i=i Pis i) 1 ^ 2 - 

We investigate three choices of detection statistic: (i) the combined matched filter 
SNR of coincident candidates before signal-based vetoes are applied (Pc lst ), (ii) the 
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Bank mass range 


2M Q < M < 35M Q 


2M Q <M< 100M q 


Statistic 


Statistic 


Found 


Statistic 


Found 




Threshold 


Injections 


Threshold 


Injections 


„rirst 
rc 


9.18 


73 


9.8 


91 


-second 
rc 


9.18 


69 


9.8 


93 


-second 
Pctt 


10.05 


27 


10.05 


85 



Table 6. Number of injections found as determined by the Neyman- 
Pearson criteria for different choices of detection statistic A and threshold A* . 
The mass range of the template bank is shown in the first row, all other parameters 
of the search as the same as those described in column (1) of Table 5. 

combined matched filter signal-to-noise ratio after the x 2 signal-based veto has been 
applied applied to coincidences (p^ ocond ), (iii) the combined effective SNR (plff° nA )- 
This statistic is only available after the second coincidence stage, since it is a function 
of matched filter SNR and the \ 2 statistic for a candidate. To set a threshold for each 
statistic we choose the highest value of that statistic NINJA data containing only noise. 
To do this, we discard all triggers within 5 s of an injected signal; the remaining triggers 
will be due to the simulated noise alone (we note that this approach is not possible 
in real data where the locations of the signals are unknown). This crude method of 
background estimation should provide us with consistent criteria for elimination of 
spurious detections. Therefore, we mark an injection as found only if it resulted in a 
trigger with statistic higher then any background trigger found in the data. 

Table 6 shows the threshold and the number of triggers found for each choice 
of statistic. It is interesting to compare the results for the low-mass search when we 
threshold on p^ ccond , rather than using a 50 ms time window to determine detected 
signals. When using the time-window method, the number of injections found by the 
low-mass search is 51, but this increase to 69 when using the the threshold method. 
Since all the injected signals lie outside the boundary of the low-mass bank, the 
coalescence time of the signals will be poorly estimated. This will result in triggers 
outside the 50 ms window, which are nevertheless arc loud enough to lie above the 
background. 

Signal-based vetoes are applied at the second stage of the inspiral pipeline and 
are used to compute p e ff ■ By comparing the number of triggers found before and after 
signal-based vetoes are applied, we can evaluate their effect on the sensitivity of the 
search. Note that we observe the same threshold for both p^ st and / o^ econd . However, 
the number of detected signals in the low-mass search is reduced by 4 as the the x 2 
veto has removed triggers where the templates are not a good match for the signals. 
More intriguing is a slight increase in the number of detected signals after the x 2 veto 
in the bank with the extended mass range (from 91 to 93). Additional investigations 
revealed that, despite having fewer triggers in each detector after the x 2 test has been 
applied, the total number of coincident triggers actually increases. This is due to the 
fact that the signal-based vetoes cause the time of the signal to be measured more 
accurately in the detectors; more triggers therefore survive the coincidence test. We 
do not observe this in the case of the low mass search. 

Finally, we turn our attention to the effective SNR statistic, defined in 
equation (8). Since the NINJA detector noise is stationary and Gaussian, the expected 
value of the x 2 1S one P er degree of freedom. Therefore, we do not expect that the 
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effective SNR will be useful in reducing the significance of loud background triggers. 
This is borne out by the fact that the statistic threshold actually increases slightly 
when using effective SNR. For the low mass search the number of signals found by 
thresholding on p c g is significantly less than when using the combined SNR statistic. 
This is to be expected as the simulated signals do not match well with the templates. 
Although the low mass templates produce candidates, these will have large values 
of x 2 since signal and template do not match well. Thus, the effective SNR will be 
smaller than the original SNR and fewer signals will be recovered above the threshold. 
This effect is less significant for the second search with a larger mass range as the 
templates provide a better match to the simulated signals. Since effective SNR has 
been a powerful statistic in real detector data, this highlights the need for further 
NINJA studies using data containing non-stationary noise transients. 

4-1-3. Inspiral- Merger- Ringdown Templates The calculation of the full binary 
binary black hole coalescence waveform accessible to ground-based detectors requires 
numerical methods. At the moment, it is not possible to accurately model a 
coalescing binary over hundreds of orbits due to the computational cost of evolutions. 
Furthermore, it is not necessary to model the entire waveform, since post-Newtonian 
gives a valid description of the system when the black holes are sufficiently separated. 
During their final orbits before merger the black holes' velocities increase and the 
post-Newtonian expansion becomes less reliable. At this stage the non-perturbative 
information contained in numerical simulations is required. A successful approach has 
been to combine analytical and numerical results to obtain full waveform templates. 
Two different families of such waveforms have been used to analyse the NINJA data: 
the effective one body (EOB) [172-175] and phenomcnological [63,67] models. 

By combining together results from post-Newtonian theory and perturbation 
theory, the EOB model [172, 173] predicts the full inspiral, merger and ringdown 
waveform. More recently, the non-spinning EOB model has been further improved 
by calibrating it to NR results, achieving high overlaps without the need to maximise 
the intrinsic mass parameters of the binary [62,64,65,68-71]. The LSC Algorithm 
Library (LAL) [149] contains two implementations of the effective one body template: 
one (called EOB) which only evolves the waveform to the light-ring frequency 

c 3 

/LR = ^73^' 

and a second (called EOBNR) which implements the full EOB waveform described 
in [65]. This template which was constructed to match the NASA-Goddard binary 
black hole simulations with mass ratios m\:mi = 1:1, 3:2, 2:1 and 4:1, however 
LAL waveforms do not yet implement higher harmonics of the signal. Both of these 
implementations were used to search for black hole binary signals in NINJA data. 

Another approach for constructing the full waveform is to "stitch" together the 
results of post-Newtonian and numerical relativity calculations. The model presented 
in [63, 67, 176] consists of matching the post-Newtonian and numerical waveforms in 
an appropriate matching regime (where both are sufficiently accurate) to obtain a 
"hybrid" waveform. This hybrid is then fit by a phenomenological model in the 
frequency domain determined entirely by the physical parameters of the system. This 
procedure has been carried out for non-spinning black holes and a two-dimensional 
template family of waveforms that attempts to model the inspiral, merger and 
ringdown stages for non-spinning binary black holes has been obtained. Each waveform 
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Template 


EOB 


EOBNR 


Phenom 


Freq. Cutoff 


Light ring 


Full waveform 


Full waveform 


Filter Start Freq. 


40 Hz 


30 Hz 


30 Hz 


Component Mass M 


10-60 


15-160 


20-80 


Total Mass M© 


20-90 


30-200 


40-160 


Minimal Match 


0.97 


0.99 


0.99 


Found Single (HI, H2, LI, VI) 


91, 64, 82, - 


97, 68, 92, 102 


92, 61, 87, - 


Found Coincidence (LIGO, LV) 


83, - 


88, 106 


81, - 


Found Second Coincidence (LIGO, LV) 


80, - 


85, 102 


80, - 



Table 7. Results of the search for NINJA signals using IMR template 
banks. There were 126 injections performed into the analysed data. The signal- 
based vetoes have little influence in the rejection of triggers, confirming their 
efficiency in separating inspiral-likc signals from other kind of glitches. 



is parametrised by the physical parameters of the system, i.e., the masses mi and m,2 
of the black holes. 

Since the EOBNR and phenomenological models provide complete waveforms, 
the search was performed to higher masses (200M Q and 16OM respectively) than 
for inspiral only searches. In principle, the search could be extended to even higher 
masses, but technical issues with the current waveform generation procedures prevent 
this. The minimum component mass was also increased, in an effort to reduce the size 
of the template bank by limiting the number of highly asymmetric signals. Finally, the 
template bank for all these searches was constructed using the standard second order 
post-Newtonian metric, and hexagonal placement algorithm [155]. At high masses, 
the parameter space metric for the full waveforms will differ significantly from the 
standard second-order post-Newtonian metric. However, the current template bank 
placement suffices for detection purposes, although probably not for good parameter 
estimation. 

The parameters of the NINJA analyses using the EOB, EOBNR and 
phenomenological waveforms are also given in Table 7. Again, the primary result is 
the number of gravitational-wave candidates found to be coincident with an injected 
signal. For the EOB model truncated at light ring, the parameters were chosen to 
match the TaylorF2 analyses described in Section 4.1.1. Therefore, it is unsurprising 
that the results are very similar to the TaylorF2 search extended to ERD (the fourth 
column of Table 5). Further details of the EOB search, and a comparison to TaylorF2 
results are available in The EOBNR results show some improvement for detecting the 
numerical relativity signals over the usual post-Newtonian or EOB waveforms. For 
the phenomenological waveforms, time windows of 120ms in single detector and 80ms 
in coincidence have been used to associate triggers to injections. These parameters 
differ from those employed in other searches to compensate for a relatively large 
observed error in the estimation of the coalescence time. By comparing the results 
with the standard post-Newtonian analyses presented in Section 4.1.1, we conclude 
that in the present case the phenomenological waveforms [63, 67] do not seem to 
provide a clear benefit over the usual post-Newtonian waveforms extended to higher 
cutoff frequency and/or to unphysical regions of the parameter space [64,72]. For 
an extended description of the search with phenomenological waveforms see [177]. 
In all cases, the signal-based vetoes have little influence in the rejection of triggers, 
confirming their efficiency in separating inspiral-like signals from other kind of glitches. 

Plots of found and missed injections for the searches are shown in Figure 9. For 
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(a) Search with EOBNR templates. 



(b) Search with phenomenological templates. 



Figure 9. Found and missed injections for the EOBNR and 
Phenomenological templates. The figures shows found and missed injections 
as a function of the injected effective distance in Hanford and the total mass. 
Left: Results for the EOBNR search. Right: Results for the search with 
phenomenological waveforms. The vertical bars mark the limits of the template 
bank used in the search. 



the most part, simulated signals in the mass range covered by the template banks are 
well recovered. Some of the missed signals at lower distance correspond to waveforms 
from simulations of spinning black holes. Since all searches make use of non-spinning 
waveforms this drop is expected. Finally, we turn to parameter estimation. Figures 10 
and 11 show the parameter recovery accuracies for the EOBNR and phenomenological 
searches respectively. In both cases, the accuracy of recovering the total mass of the 
simulations is greatly improved over TaylorF2 waveforms shown in Figure 8. This is 
likely related to the increased mass range of the searches, as well as the use of full 
waveforms. The timing accuracy for EOBNR is comparable with the TaylorF2 results, 
while for the phenomenological waveforms, the known timing bias affects the results. 

Both the EOBNR and phenomenological models will be improved in the future. 
Further accurate EOBNR models have already appeared in the literature [65,68- 
71] since the time the EOBNR model used in this analysis was implemented, 
and extensions to include spin and eccentricity are under development. There 
are a number of obvious improvements in the phenomenological waveforms that 
can be made: Calculating the parameter space metric for the phenomenological 
waveforms would enable the use of an optimal template bank and allow for improved 
coincidence algorithms. The construction of the phenomenological waveform model 
can itself be significantly improved by extending the fitting to higher mass ratios 
and spins, quantifying the error on the phenomenological parameters, matching to 
post-Newtonian theory as early as possible and including higher order modes in the 
waveform. The results of the NINJA analysis also demonstrate a clear need to improve 
accuracy in measuring the end time of the signal. This is not straightforward, however, 
since there is no clear definition of the time of merger for the phenomenological 
waveforms or the numerical signals [75]. Work on the improvements to both the 
EOBNR and phenomenological searches are being made, and will be applied in and 
guided by future NINJA projects. 
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Figure 10. Parameter accuracy for EOBNR templates. Left: Accuracy 
with which the total mass is recovered. The template bank covers the region 
30M Q < M < 200M Q , hence the mass of injections with M > 200 Mq are always 
underestimated. Most of the injections with total mass less than 200 Mq were 
recovered with a mass accurate within to a few tens of percent , demonstrating that 
the EOBNR templates are more faithful to the injected signal than the TaylorF2 
templates shown in Figure 8. Higher mass injections are necessarily recovered 
with underestimated total mass, because the template bank did not cover the 
entire simulation region. Right: Accuracy of determining the coalescence time of 
the injections. The end time for injections with total mass less than 200 Mq was 
typically recovered to within a few milliseconds. The end time for injections with 
total mass above 200 Mq (outside the range of the template bank) was typically 
recovered to within 10 or 20 milliseconds. 
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Figure 11. Parameter accuracy for phenomenological templates. Left: 

Accuracy with which the total mass is recovered. The total mass is typically 
recovered within 20%, for signals within the template space. For higher mass 
injections, there is an inevitable underestimation of the mass due to the limited 
reach of the template bank. Right: Accuracy of determining the coalescence 
time of the injections. The timing plot shows the systematic offset discussed in 
the text. 
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4- 1-4- Ringdown Templates As described in section 4.1, ringdown templates can be 
computed using black hole perturbation theory and so matched filtering can be used 
to search for these signals. Ringdown templates are exponentially damped sinusoids 
parametrised by the ringdown frequency / and quality factor Q. The LSC ringdown 
search pipeline [160] has been used to filter the NINJA data against a bank of ringdown 
templates with frequencies between 50 Hz and 2 kHz, and quality factors between 2 
and 20. The bank had a maximum mismatch of 3% and contained 583 templates. A 
lower-frequency cutoff of 45 Hz was applied when filtering the NINJA data generated 
with the LIGO noise curves and 35 Hz for data with the Virgo noise curve. The 
goals of these analyses were to ascertain the detectability of the injected numerical 
waveforms using ringdown templates at single and coincident detector levels and the 
accuracy with which the final black hole parameters can be estimated. The current 
searches use single-mode templates. The waveforms described in this paper are known 
to contain higher order multipoles. The potential effects of ignoring these in the search 
are discussed in Ref. [178] (see in particular Fig. 8 in there). 

An injection is defined as found if a set of coincident triggers lies within 10 ms 
of the peak time of the injection (as specified in the contributed numerical data). If 
more than one set of coincident triggers satisfies this criterion, that with the largest 
value of J2i Pi i s selected, where p t is the signal to noise ratio in the i tn detector. 
Of the 126 injections made into the three simulated LIGO detectors, 45 were found 
in triple coincidence, 24 in HI and LI (only), and 7 in HI and H2 (only). Figure 12 
shows the distribution of found and missed injections for this analysis. The ringdown 
frequency and quality is computed via the Echeverria formulae [179]: 



More recent and accurate fits for a variety of modes are listed in the Appendices 
of Ref. [180]. The final black hole mass M and spin a can be computed from the 
component masses and spins of the numerical simulation, as described in [65] and 
[181], respectively. See also Refs. [51, 182] for a discussion and comparison of different 
numerical techniques to perform the necessary fits. 

As expected, we see that in general, the closest injections (measured by effective 
distance -D c ft), defined in equation 10) were found in triple coincidence, those with a 
large Livingston effective distance were found in HI and H2 only, while those with a 
large Hanford effective distance were not found in H2, and the furthest injections were 
missed in at least two detectors. The plots show that there are three missed injections 
which, given their frequencies and effective distances, we would have expected to find. 
However, all three of these arc (non-spinning) injections with mass ratio of 4:1, and 
thus the energy emitted in the ringdown is less than would be emitted by a binary of 
the same total mass but with a mass ratio of 1 [51]. This is not taken into account in 
the calculation of effective distance. 

Equations (12) and (13) can be inverted to calculate M and a from the template 
parameters / and Q of a given gravitational-wave candidate. Figure 13 shows the 
accuracy with which the ringdown search measures the mass and peak time of the 
injected signals. Given that mass is radiated during the ringdown phase (the exact 
amount depends on the initial mass ratio) one would expect the measured mass to 
underestimate the mass of the signal, and hence the data points would lie below the 
diagonal. However, the recovered frequency is systematically underestimated due to 
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Figure 12. Distribution of injections found and missed by the ringdown 
pipeline. The left figure shows the effective distance of the injected signal in the 
LIGO Hanford Observatory as a function of the predicted ringdown frequency. 
The right figure shows the effective distance of the injected signal in the LIGO 
Livingston Observatory as a function of the total initial mass of the signal. 
The figures show signals found in triple coincidence (blue crosses), in double 
coincidence in H1H2 (green stars), in double coincidence in H1L1 (cyan stars), 
and missed (red circles). 
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Figure 13. Accuracy of measuring the ringdown parameters. The left 
figure shows the detected ringdown mass versus total injected mass for all found 
injections. The right figure shows the difference between the time of injected 
waveform peak amplitude and the start time of the ringdown as found by the 
search. 



the presence of the preceding inspiral, leading to an overestimation of the mass. The 
peak time of the signal is measured with similar accuracies to the coalescence time 
measured by the TaylorF2 templates described in Section 4.1.1. The three data points 
with a large time difference and masses lying in the range 80 and 110 M Q are part 
of the PU_T52W non-spinning, equal mass group where the peak amplitude occurred 
early in the waveform (i.e prior to the merger). 
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4-2. Search pipelines to detect un-modelled waveforms 

Several algorithms exist to detect gravitational wave transients with minimal 
assumptions on their origin and waveform; these techniques are often referred to as 
burst searches. Burst searches do not use templates and instead target excesses of 
power in the time-frequency plane. The LSC and Virgo collaborations have developed 
several burst search algorithms which use different transformations for the generation 
of time-frequency data maps. The identification of coherent signatures across multiple 
detectors has proven to be very effective at suppressing false alarms. 

Since they do not assume a template and they target short transients burst 
searches are suited for the detection of the merger phase of the coalescence. They 
have the potential to probe a large parameter space, inclusive of spin and ellipticity, 
at no additional computational cost. For this reason, the NINJA data was analysed 
by two burst algorithms: Q-pipeline and HHT. 

The Q-pipeline [183, 184] is one of the algorithms used in to search for burst 
sources in LIGO's fifth science run [185]. It is a multi-resolution time-frequency 
search for statistically significant excess signal energy, equivalent to a templated 
matched filter search for sinusoidal Gaussians in whitened data. The template bank 
is constructed to cover a finite region in central time, central frequency, and quality 
factor such that the mismatch between any sinusoidal Gaussian in this signal space 
and the nearest basis function does not exceed a maximum mismatch of 20% in energy. 
For the purpose of the NINJA analysis, and to explore detectability and parameter 
estimation, the Q-pipcline analysis was focused on the detection efficiency at the single 
detector, for all four detectors, using a nominal SNR threshold comparable to the one 
used in the matched filter searches. 

The Hilbert-Huang Transform (HHT) [186, 187] is an adaptive algorithm that 
decomposes the data into Intrinsic Mode Functions (IMFs), each representing a unique 
locally monochromatic frequency scale of the data. The original data is recovered by 
constructing a sum over all IMFs. The Hilbert transform as applied to each IMF 
unveils instantaneous frequencies and amplitudes as a function of time, thus providing 
high time-frequency resolution to detected signals without the usual timc-frcquency- 
uncertainty as found in basis set methods like the Fourier transform. 

In this section we briefly describe how the algorithms were applied and highlight 
their performance, while Section 4.3 compares the performance of the burst searches 
to the matched filtering algorithms. 

4-2.1. Q-pipeline The simulated LIGO and Virgo data streams were filtered by the 
Q-pipeline [210] with the same configuration used in the LSC S5 burst analysis [185]. 
Data is processed in 64 s analysis blocks with frequency range 48-2048 Hz and Q range 
3.3-100. The resulting triggers, once clustered, indicate a time-frequency interval and 
a significance of the excess power in that time-frequency tile. This significance can be 
easily converted into the signal-to-noise ratio of a matched filter with sine-Gaussian 
templates. 

Figure 14 shows the distribution of found and missed injections in the HI 
detector as a function of the total mass and the matched filter SNR of the injected 
waveforms for Q-pipcline (left) and the EOBNR matched filtering search described in 
section 4.1.3, where both detection thresholds are set at SNR = 5.5. At the single 
detector level, the two algorithms have comparable performance. Figure 15 shows the 
the central frequency of the most significant tile reported by Q-pipeline versus the 
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Figure 14. Distribution of injections found by the Q-pipeline (left) in 
the LIGO Hanford 4 km detector, to be compared to the distribution 
for the EOBNR matched filtering search. The total mass of the injected 
signal is shown on th x-axis and the optimal matched filer signal-to-noise ratio 
of the injection is shown on the j/-axis. Circles show found injections and crosses 
show missed injections. At the single-detector level, with the same SNR threshold, 
Q-pipelinc and the EOBNR search have comparable performances. 




Figure 15. Parameter measurement of signals using the Q-pipeline. The 

left plot shows the measured frequency of the gravitational- wave candidate versus 
the injected ringdown frequency / r i ng computed by equation (12). The right plot 
shows the difference between Q-pipcline measured frequency and / r i ng , against the 
total mass of the injection. These plots show the algorithm preferentially triggers 
on the portion of the signal that is in the most sensitive band of the detector 
(50-200 Hz). This is consistent with the behavior of the ringdown search; see 
figure 18 for an event-by-event comparison. 



ringdown frequency computed from equation (12) and the difference between these 
two frequencies as a function of mass. These results demonstrate that the Q-pipelinc 
preferentially detects the portion of the signal that is in the most sensitive frequency 
band of the detector (50-200Hz): the ringdown for higher masses, or the inspiral for 
lower masses. 

4-2.2. Hilbert-Huang Transform An automatic, two-stage HHT pipeline was applied 
to the NINJA data to detect and characterise the injected signals. Since the HHT 
pipeline [188] is a new development its application to NINJA data was restricted to 
the simulated 4 km LIGO detectors HI and LI. The data was pre- whitened and a 
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Figure 16. Distribution of found and missed injections for the Hilbert— 
Huang pipeline. The left figure shows the results of the search on the simulated 
LIGO Hanford 4 km detector, and the right figure shows the result for the 
simulated LIGO Livingston 4 km detector. Detected signals are shown in black 
and missed signals in white, as functions of the injected matched filter signal-to- 
noise ratio and total mass. 



1000 Hz low pass zero-phase finite impulse response filter was applied prior to use 
of the HHT. In the first detection stage, the instantaneous amplitudes from each 
detector in turn are divided into blocks with similar statistical properties according 
to the Bayesian Block algorithm [189]. These blocks are then scanned for excess 
power, with triggered blocks yielding start and end times, thus coincidences between 
detectors, the maximum frequency, and the signal-to-noise ratio (SNR) of the signal. 
The second characterisation stage computes the instantaneous frequencies, detailed 
time-frequency-power maps and time-frequency-power cluster-enhanced maps for the 
region of data containing the signal identified in the first stage. The overlap of the 
individual cluster-enhanced maps is used to estimate the time lag between the signals 
in the detectors and to construct a coherent addition of the two detector data streams 
used in a final characterisation of the signal. 

The excellent resolution of the HHT in time and frequency was used to reject false 
events due to overly short triggers, failed coincidences or mismatched time- frequency 
ranges [188]. We identified the latter as a powerful veto tool which could be used 
to improve the sensitivity of gravitational-wave data analysis pipelines. We were 
ultimately able to identify 80 events in coincidence, as shown in Figure 16. Three 
of these candidate events were determined to be false positives when the candidate 
events were compared to the list of injected signals. Out of the 50 missed events, 39 
have injected SNR < 10, five have injected SNR < 10 in one detector and SNR > 10 
in the other, and six had injected SNR > 10 in both detectors. We therefore reason 
that most of the of missed events are low SNR cases in which no blocks were triggered. 
Most of the cluster-enhanced maps show the time-frequency evolution of the signal 
with high precision (see Figure 17 for one particular example). Time lag estimates 
and coherent additions show strong potential and can be seen as proof of principle, 
but need further refinement to work reliably in an automatic pipeline. SNR estimates 
are difficult since only the burst region or diverse fragments of the signal are visible 
in our search. We refer to [190] for further details of this analysis. 
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Figure 17. The cluster-enhanced time-frequency maps of the BAM 
HHB S00 signal injected with a total mass of 69.8A/q. HI is shown in the 
left panel and LI in the right panel. We clearly see the burst part of the signal, 
thus the actual merger. The ringdown and inspiral may be obscured by noise at 
low SNR. 



4-3. Comparison of Inspiral- Burst- Ringdown Results 

In this section we consider comparisons between several of the search pipelines 
described in the preceding sections. The performance of a pipeline depends on 
many parameters, such as the signal-to-noise thresholds, the trigger coincidence and 
coherence tests, signal-based vetoes tests, allowed false alarm rates, etc. Search 
pipelines are tuned to suppress false alarms while preserving detection efficiency 
and tuning decisions can dramatically affect the relative performance of one pipeline 
versus another. Given the limited scope of the NINJA data set, a comprehensive 
comparison of pipelines was not possible. However, to make a first-order comparison 
between search pipelines applied to the NINJA data, disentangled from pipeline tuning 
decisions, we compared the number of injections found in a single detector at a fixed 
matched-filter signal-to-noise threshold. 

Table 8 reports the number of injections found in single interferometers and in 
multi-interferometer networks, using triggers from the Q-pipclinc burst algorithm 
(which targets the merger by match filtering to sine-Gaussian templates), matched 
filter to ringdown templates, matched filter to inspiral templates and matched filter 
to non-spinning, full coalescence EOBNR waveforms. For all algorithms the same 
nominal threshold of SNR > 5.5 was imposed. In addition, the number of detected 
injections in AND (injections detected simultaneously from multiple algorithms) and 
OR (injections detected by at least one algorithm) in an inspiral-mcrgcr-ringdown 
analysis is reported. 

The statistics of this sample is too small to make inferences on which pipeline 
performs better in which parameter region; a more systematic study is needed. 
Furthermore, since the NINJA data contains only Gaussian noise and signals and 
so this comparison does not take into account noise transients which may cause false 
detections. Nevertheless, we have an indication that all pipelines have comparable 
chances to find these injections. The differences between pipelines are in the accuracy 
with which they can measure the parameters of the signal. 

Figure 18 compares the accuracy with the Q-pipclinc and ringdown searches 
measure the frequency of the signal and Figures 19 compares the accuracy with these 
pipelines measure the total mass of the injected signal. Both the ringdown and Q- 
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Table 8. Number of injections found with SNR > 5.5. This table 
takes into account for each detector only signals with injected SNR > 5.5. 
If the same injection is found in more than one detector or algorithm, it is 
counted as coincidence. The inspiral triggers are from the 2 PN TaylorF2 
templates, terminated at /isco (total mass in 20 — 90 Mq) and from the EOBNR 
waveforms. Note that, in an actual search, different thresholds may be needed 
for different pipelines, depending on the false alarm rate and the morphology of 
noise transients, so this is strictly a comparison between search techniques on the 
limited NINJA data set, not of the performance of full pipelines on actual data. 




Figure 18. Comparison of the frequencies measured by the Q-pipeline 
and ringdown searches. The x-axis shows the difference between the central 
frequency reported by the Q-pipeline and the frequency of the injected ringdown 
/ring- The y-axis shows the difference between the frequency measured by the 
ringdown search and the injected ringdown frequency. 
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Figure 19. Comparison of the injected and recovered masses for the 
Q-pipeline and ringdown searches. For both searches we start from the 
frequency reported by the algorithm and compute the corresponding mass, using 
equation (14), assuming that the measured frequency is /iscOi /lr, and /erd- 
Both algorithms measure a frequency somewhere between the light-ring and 
ringdown frequencies. 




Figure 20. Comparison of the injected and recovered masses for the 
TaylorF2 and EOBNR searches. The left plot shows the comparison for 2.0 
order post-Newtonian TaylorF2 templates terminated at /isco ana - the right plot 
shows the comparison for EOBNR templates. 



pipeline searches report a frequency not a mass, so to compare to the injected and 
detected total masses we must calculate a mass from the frequency of the candidate. 
There is no unique way to do this, and but we can calculate the total mass under 
the assumption the algorithm is detecting a given portion of the waveform. Figure 19 
shows the result of using the formulae for ISCO, light-ring and effective ringdown 
frequencies to compute the total mass [64, 173]: 

c 3 c 3 0.5967c 3 nA . 

/lSC °-tV6V^? /LR -373V^? hw '^GM- (U) 
Both burst and ringdown code preferentially detect the portion of the signal that is in 
the detectors' most sensitive band. Figure 19 shows that for both algorithms, at the 
lowest injected masses this correspond to fisco, but as the injected mass increases, 
the algorithms trigger between light ring and ringdown, as expected. 

This comparison is more straightforward for matched filtering codes that use a 



Results from the first NINJA project 



40 



template with a given mass. Figure 20 shows the detected total mass against the 
injected total mass for (i) 2.0 post-Newtonian order TaylorF2 templates terminated 
at /iscCb with templates in the mass range 20-90 M Q and (ii) EOBNR templates 
with masses in the range to 30-200 M Q . The TaylorF2 templates significantly 
underestimate the masses of the injected signals, due to the fact that most of the 
injected signals lie outside the template bank. The EOBNR search, with its larger mass 
range and inspiral-merger-ringdown waveforms, measures the masses of the injected 
signals with better accuracy. 



4- 4- Bayesian pipelines 

Bayesian inference [191] is a powerful means of extracting information from 
observational data based on the calculation of posterior probabilities and probability 
density functions. Computation of these quantities is expensive and so these 
algorithms are not typically used to search for candidate events in detector data. They 
are, however, useful in the closer study of candidates identified by the search pipelines 
described in Sections 4.1 and 4.2. This section describes the results of applying two 
different Bayesian inference algorithm to the NINJA data. The first is designed to 
estimate the parameters of the signal assuming a gravitational wave is present in the 
data. The second calculates the confidence in the presence of the signal, quantified by 
the odds ratio between the signal and noise models of the data. 

Both approaches require the calculation of the posterior probability-density 
function (PDF) on the parameter space of the signal, given the data d, which is 
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in the presence of Gaussian noise with power spectral density S n (f), where p(9) is the 
prior probability density of the parameters 9 and h(9) is the model used to describe 
the signal [192]. 

A Markov-Chain Monte-Carlo approach [193] was used to coherently analyse data 
from multiple detectors in order to evaluate the posterior PDFs. This technique 
stochastically samples the parameter space in a search for the parameters that best 
match the observed data; it does so by attempting a random jump from the current 
set of parameter values to a new one, then deciding whether the jump should be taken 
by comparing the likelihoods of the old and new locations in parameter space. In this 
way, one simultaneously searches for the set of parameters that yield the best fit to 
the data, and determines the accuracy of the parameter estimation. 

Bayesian model selection, based on a different Monte-Carlo technique known 
as nested sampling [194], was employed as a tool to measure the confidence of a 
detection using different waveform families. This approach requires the calculation 
of the marginal likelihood of the signal and noise models, obtained by computing the 
integral of the posterior PDF J p{9)p(d\9)d9 to find the total probability of the model. 
It was possible to calculate this integral for the nine-parameter model of a non-spinning 
binary coalescence signal described coherently in multiple interferometers. The ratio 
of likelihoods of the signal and noise models is known as the "Bayes factor," and is used 
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to multiply prior odds, giving the posterior odds ratio between the two models, taking 
into account the observational data; in turn, the value of this Bayes factor corresponds 
to the level of confidence in the detection. As a straightforward by-product of the 
nested-sampling algorithm, it is also possible to infer the maximum-likelihood values 
of the parameters of the detected signals; this was used to obtain further information 
on the ability of different waveform approximants to recover the source parameters. 

4-4-1- Parameter Estimation Using Markov Chain Monte Carlo A selection of 
injected numerical signals were analysed with a Markov-Chain Monte Carlo (MCMC) 
code [195, 196]. The signal model was based on waveforms with phase evolution at 1.5 
post-Newtonian order and leading-order amplitude evolution. Parameter estimation 
was successful on NINJA injections with relatively low total mass in which the inspiral 
contained a significant fraction of the total signal-to-noise ratio. For high-mass 
injections, the algorithm attempted to match the merger and ringdown portions of 
the injected signal to inspiral templates, resulting in poor parameter estimation. 

The post- Newtonian waveforms used in this analysis include the spin of the larger 
body mi, allowing us to use the analytical simple-precession waveform [197]. The 
parameter space thus consists of twelve independent parameters: 

9 = {M, ?y, a, 8, ip, l, d, a spin , k, 4> c , a c ,t c }, (16) 

where M. and 77 are the chirp mass and symmetric mass ratio, respectively; a (right 
ascension) and S (declination) identify the source position in the sky; the angles ip and l 
identify the direction of the total angular momentum of the binary; d is the luminosity 
distance to the source; < a sp i n = Si/m\ < 1 is the dimensionless spin magnitude; k 
is the cosine of the angle between the spin and the orbital angular momentum; and 
4> c and a c are integration phase constants that specify the gravitational-wave phase 
and the location of the spin vector on the precession cone, respectively, at the time of 
coalescence t c . 

The MCMC algorithm used for the NINJA analysis was optimised by including a 
variety of features to efficiently sample the parameter space, such as parallel tempering 
[196]. This MCMC implementation can be run on a data set from a single detector, or 
on data sets from multiple detectors. In the latter coherent search among all 

detectors significantly improves the determination of source position and orientation 
[195,196]. 

The MCMC code was run on a selection of injected signals in the NINJA data. It 
was found that although the MCMC runs are clearly able to detect a signal whenever 
the inspiral contains a sufficient signal-to-noise ratio (SNR), they were unable to 
correctly determine the signal parameters for many injections. For the high masses 
typical of most NINJA injections, the SNR is dominated by the merger and ringdown, 
so that the inspiral-only templates tried to match the merger and ringdown portions of 
the injected waveform. Typically, it is found that in such cases the time of coalescence 
is overestimated since the injected ringdown is matched to an inspiral; the chirp mass 
is underestimated since the merger/ringdown frequency is higher than the inspiral 
frequency for a given mass, so that matching them to an inspiral requires the mass 
to be lower; the mass ratio is underestimated, which allows the waveform to contain 
more energy in the narrow frequency band corresponding to quasi-normal ringing; 
and the spin rails against the upper prior of 1 since the innermost stable circular 
orbit frequency is highest for an inspiral into a maximally-spinning Kerr black hole. 
We tried to circumvent the problem of matching to the merger and ringdown by 
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introducing more-restrictive priors on spin and/or rj. These efforts stili failed when 
the total masses were too high, but were successful in the case of lower total masses 
and longer inspiral signals. 
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Figure 21. Results of the MCMC analysis for one of the injected 
NINJA waveforms. The top row shows the marginalised PDFs for M, r), 
a spim an d d produced by a two-detector MCMC analysis on an injected non- 
spinning equal-mass SpEC Cornell/Caltech waveform (true values Ai = 15.6Mq, 
rj = 0.25, t c = 4.7223 s, and a spin = 0). Middle row: the same PDFs (but with 
Gspin replaced by t c ) for a three-detector run with constrained spin on the same 
injection. Bottom row: two-dimensional PDFs for the sky position with the 2- 
detector run on the left and the three-detector MCMC run on the right; the 1-a, 
1-a and 3-cr probability areas are displayed in different colours, as indicated in 
the top of each panel. Dashed lines denote the true values of injected parameters. 

Figure 21 shows the PDFs produced by runs on an injected equal-mass non- 
spinning SpEC Cornell/Caltech waveform with M = 15.6 M@. This particular injection 
was chosen because it had the lowest total mass, and SpEC waveforms typically have 
more inspiral cycles; runs on other injections show comparable results, with the general 
trend that the higher the total mass (and, thus, the lower the relative fraction of the 
SNR contributed by the inspiral), the poorer the parameter estimation becomes. 

Data from two detectors, HI and LI, were used to compute the PDFs shown in 
the top row of Figure 21. We used wide, flat priors for intrinsic parameters (e.g., 
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M G [2M Q ,100M Q ], 77 S [0.03,0.25], a spin e [0,1]). We find that the values of the 
intrinsic parameters are not determined very accurately. In particular, the spin rails 
against the upper bound of 1 while r\ is underestimated, as expected. We find that the 
sky location is nevertheless constrained to an arc of a ring containing the true value; 
the 2-cr (~ 95%) sky area of the ring shown in the bottom left of Fig. 21 is ~ 10000 
square degrees. 

In the middle row of Fig. 21, we plot the PDFs based on data from three detectors: 
HI, LI, and VI. The spin parameter was constrained to its true value a sp i n = 
for this run. This had the effect of forcing the MCMC to match the inspiral only, 
significantly improving the resolution of other parameters: for example, the PDF of 
r\ now rails against 0.25, which is its true value. The chirp mass is still somewhat 
underestimated; a higher SNR may be necessary to improve the mass determination. 
Promisingly, it was found that the sky location is constrained to a smaller patch 
on the sky: the 2-cr sky area in the bottom right of Fig. 21 is 6300 square degrees. 
In fact, the sky localisation is even better when the spin parameter is allowed to 
vary, allowing the MCMC to use the SNR contributed by the ringdown; removing 
the spin-parameter constraint reduces the 2-cr sky area to 2750 square degrees. This 
ability to determine the source position will be significant in any future searches for 
electromagnetic counterparts of gravitational-wave triggers. 

We hope that in the future it will be possible to test MCMC codes on numerical 
signals in a lower mass range, where the inspiral portion would dominate the SNR, 
so that inspiral-only templates are not at a significant disadvantage. Meanwhile, we 
have recently implemented templates at 3.5 PN order in phase that include the spin 
for both members of the binary, thus improving the accuracy of parameter estimation 
and increasing the range of applicability of our code. 



Bayesian Model Selection Pipeline The primary goal of this analysis was 
to investigate the performance of different template families on the confidence of 
detection of the injections contained in the NINJA data set. The approach to this 
problem used a method described in [198, 199], which can be summarised as follows. 
Two hypotheses arc under consideration: (i) Hn — The data {dk} are described by 
(Gaussian and stationary) noise only: dk = and (ii) Hs — The data {dk} are 
described by (Gaussian and stationary) noise {n&} and a gravitational wave signal 
{h k a \9)}, according to a given approximant a, where 9 represents the vector of the 



(unknown) signal parameters: dk = n k + h k a \9). 
calculated by performing the integral 



The marginal likelihood of Hs is 



p({d~k}\H s ,a) 



p(9)p({d k }\Hs,a 1 9)d9. 



(17) 



The ratio of probabilities or "odds ratio" of the two models is 

P(H s \{d k },a) 



O 



SN,a 



P(H N \{d k }) 
P(Hs) 



P(H N ) 



P({d k }\H s ,a) 
P({d k }\H N ) 



= VB 



SN> 



(18) 



where V is the prior odds ratio and B s a ^ is the Bayes factor. 
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In this analysis the model includes the response of all available simulated detectors 
(LI, HI, H2 and VI) coherently, and the gravitational waveforms are calculated using 
function in the LAL library [149]. For the gravitational waveform, two different 
approximants were considered: the standard (2 PN) TaylorF2 waveform family, 
with inspiral truncated at /isco, an d the phenomenological inspiral-mcrgcr-ringdown 
IMRPhcnA templates described in [63]. In each case the waveforms were truncated 
at low frequency at 30 Hz. 
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Figure 22. The distribution of Bayes factors obtained by running the 
algorithm with TaylorF2 approximants on signal-free data sets. The 

pipeline was run coherently on the four detectors, and the same was done for the 
results reported in this figure. The total number of trials is 2000 and the highest 
value of the Bayes factor is log 10 B = 2.77. 



The choice of priors in the analysis was as follows. For the TaylorF2 approximant, 
the prior on 9 was uniform on chirp mass, symmetric mass ratio and distance, within 
the following limits: time of coalescence in a window ±0.5 s around the actual t c of 
each injection, symmetric mass ratio rj in the range 0.01 < r/ < 0.25, chirp mass 
within the bounds set by r\ and the total mass in the range 50A/q < M < 150Mq, 
and distance 1 < D < 500 Mpc. The parameters for orientation and position on 
the sky of the binary were allowed to vary over their entire angular ranges. For 
the IMRPhenA approximant, the limits were identical on all the parameters, with 
the exception of the total mass whose upper boundary was set to 475M Q . We also 
calculated the Bayes factors for data segments containing no signal in order to estimate 
the background distribution of Bayes factors. Figure 22 shows the distribution of 
Bayes factors (using TaylorF2 approximants) when running the analysis algorithm on 
portions of data without injections: 2000 trials were carried out, with a maximum 
value of log 10 B = 2.77. If interpreted as a threshold value on the Bayes factors to 
decide whether a signal has been detected or not, it corresponds to a false alarm of 
0.05%. The distribution obtained with the IMRPhcnA approximant is very similar. In 
the analysis, a range of "detection thresholds" on log 10 B was considered, in order to 
explore how the two different approximants (and the algorithm) respond to different 
numerical relativity injections. 

Figure 23 and Table 9 summarise the main results. The left panel of Figure 23 
shows the value of the Bayes factor computed for the two approximants as a 
function of the coherent four-detector signal-to-noise ratio at which the waveforms 
were injected. For all the injections, IMRPhenA approximants return a Bayes 
factor which is (significantly) larger than TaylorF2 approximants. This is not 
surprising, as the TaylorF2 waveforms do not contain the merger and ring-down 
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Figure 23. Comparison of the Bayes factors for TaylorF2 and 
IMRPhenA approximants. Left: The values of the Bayes factors obtained in 
the analysis of the NINJA data set as a function of the optimal coherent (Ll-Hl- 
H2-V1) signal-to-noise ratio at which the signals were injected into instrument 
noise. The solid (black) dots and the open (red) circles refer to the Bayes factors 
obtained by using the IMRPhenA and TaylorF2 approximants, respectively. 
Right: The cumulative number of injections recovered as a function of the 
Bayes factor. The thin (red) solid line corresponds to the results obtained by 
using the TaylorF2 approximant, whereas the thick (black) solid line refers to the 
IMRPhenA approximant. A threshold of log 10 Bsn = 3 has been used. 



Threshold 


Number of detected injections 


(logio B) 


TaylorF2 


IMRPhenA 


3 


69 


112 


5 


61 


107 


10 


43 


104 


30 


28 


89 


100 


10 


58 



Table 9. The number of detections in the NINJA data with model- 
selection pipeline as a function of the Bayes-factor threshold using 
TaylorF2 and IMRPhenA approximants. The total number of injections 
was 126. See also Figure 23. 



portion of the coalescence and are truncated at /isco- Figure 23 (right panel) 
shows the number of injections that are recovered at a given Bayes factor (or above). 
Once more the effectiveness of IMRPhenA approximants is striking compared to 
the TaylorF2 waveform family. These results are broadly in agreement with the 
matched- filter analysis carried out with inspiral-merger-ringdown waveforms described 
in Section 4.1.3. 

The nested-sampling algorithm used for model-selection can also be used for 
parameter estimation. In particular, one can generate in a straightforward way the 
maximum likelihood estimate of the recovered parameters, and have an indication of 
the statistical errors on such values. For simplicity, in this analysis we identified the 
statistical errors (the error bars in Figures 24 and 25) with the region of parameter 
space in which the likelihood values were not lower than a factor e with respect 
to the maximum likelihood. The results for chirp mass, total mass and the two 
coordinates of the source in the sky - latitude and longitude - are shown in Figures 24 
and 25; here we restrict to only the IMRPhenA approximant and to all the signals 
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that yielded log 10 B > 3. The results for the masses show a behaviour that is 
qualitatively consistent with the results obtained using a matched-filtering analysis, 
see e.g. Figures 10, 11, and 20. The total mass is (in most of the cases) systematically 
underestimated, although for 34 injections the recovered values were consistent with 
the injected total mass. These injections correspond in all cases to waveforms with 
(near) zero eccentricity and in 21 (out of 34) instances to non-spinning waveforms. 
We have also checked that the errors on the masses do not show any significant 
correlation with the value of the Bayes factor at which the injections were recovered 
or the injected signal-to-noise ratio. However, despite the systematic errors on the 
physical parameters, the sky location is on average fairly well determined. This is 
most likely due to the fact that there is enough information in the (source-location 
dependent) time of arrival of the signals at different instrument sites to recover 
meaningful information about the position of the source in the celestial sphere. This 
is currently under careful investigation and more details about this and other aspects 
of the analysis can be found in Rcf. [200] . 




Figure 24. Comparison of the recovered mass parameters for the 
IMRPhenA approximant. Left: The recovered (maximum likelihood) values 
of the chirp mass as a function of the injected values. Right: The recovered 
(maximum likelihood) values of the total mass as a function of the injected values. 
The IMRPhenA approximant was used with a threshold of log 10 B$n = 3. 



5. Conclusion 

The NINJA project was conceived as a first step towards a long-term collaboration 
between numerical relativists and data analysts with the goal of using numerical 
waveforms to enhance searches for gravitational waves. NINJA is unique in that it 
focused on running existing gravitational-wave search algorithms on data containing 
waveforms obtained from numerical simulations. Since this constitutes the first such 
analysis, the scope of the project was deliberately kept somewhat modest: restrictions 
were placed on the number of waveforms to be submitted by each numerical group, 
no attempt was made to include transient noise sources in the data and only a limited 
number of simulated signals were produced for the data analysis. This helped to 
encourage significant involvement from both the numerical relativity and data analysis 
communities, with ten numerical relativity groups providing waveforms and data- 
analysis contributions from nine different groups. 
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longitude (rad) latitude (rad) 

Figure 25. The recovered source-location parameters for the 
IMRPhenA approximant. Left: The recovered (maximum likelihood) values 
of the source longitude as a function of the injected values. Right: The recovered 
(maximum likelihood) values of the source latitude as a function of the injected 
values. The IMRPhenA approximant was used with a threshold of log lc) B$ at = 3. 



Communication between the data analysis and numerical communities has been 
smooth and fluent during the course of the NINJA project. The format described 
in [76] provided a good starting point from which to interchange data between the 
communities. As the project was being developed, several improvements were made 
to the format, which we expect will continue evolving as more experience is gained 
with a broader family of waveforms, including those containing matter. 

The limited number of signals in the NINJA data makes it dangerous to draw 
strong conclusions from the comparison of different waveform families and different 
search methods. Overall, it is clear that many of the data analysis methods were 
capable of detecting a significant fraction of the simulated waveforms. This is 
immediately significant as several of the analyses performed are routinely used in 
searches of the LIGO and Virgo data. However, since the NINJA data set did not 
include the type of non-Gaussian transients seen in real gravitational-wave detector 
data, it is difficult to translate the efficiencies observed here into statements about 
LIGO or Virgo sensitivity. 

NINJA has demonstrated that more work is required to measure the parameters 
of signals in detector data. Parameter estimation is poor for most pipelines, and 
several methods tend to associate a candidate event to that part of the waveform 
which lies in the most sensitive band of the detector. For example, in a search with 
inspiral only templates, the ringdown of a high mass black hole which occurs at around 
100 Hz might be picked up. This will lead to poor estimation of both the binary's 
mass and coalescence time. Similarly, the un-modelled burst searches will correctly 
identify the signal but, without knowing which part of the coalescence it corresponds 
to, have difficulty providing accurate parameters. There is some evidence that using 
full inspiral-merger-ringdown waveform templates alleviates this problem, as well as 
evidence that estimation of the sky location of the signal is largely independent of the 
mismatches between simulated and template waveform. These are all issues which 
warrant further investigation. 

We hope that this work will provide a foundation for future analyses, and plans are 
envisioned to continue and extend the NINJA project. Several suggestions have been 
made to broaden this work and make it more systematic: in addition to expanding 
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the parameter space explored by numerical simulations, two crucial steps will be 
to construct hybrid analytic-numerical waveforms (which will allow a lower range 
of masses to be injected) and to consider data containing non-stationary noise. It 
would also be natural to include other waveform families, such as supcrnovae or 
binary mergers comprising one or two neutron stars. Subsequent NINJA projects 
could provide a noise-free data set for tuning parameter estimation and measurement 
pipelines and release "training" and "challenge" data sets, as has proven successful 
in the Mock LISA Data Challenges [201,202], in which the parameters of the 
waveforms are known and unknown to the analysts, respectively. The numerical 
data sets may also be useful for efforts aimed at using the best-available waveforms 
to explore and develop LISA data analysis approaches and in evaluating parameter 
estimation accuracy for LISA. These efforts, as carried out by the Mock LISA Data 
Challenge Task Force and the LISA Parameter Estimation Task Force, are summarised 
in Ref. [203,204]. 

However future analyses progress, it is clear that a significant amount remains to 
be learned from collaborations between the numerical relativity and gravitational- wave 
data analysis communities. 
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Term Meaning 

ADM Arnowitt-Deser-Misner. 

ASCII American Standard Code for Information Interchange. 

BAM Bifunctional Adaptive Mesh code developed at University of Jena. 

BBH Binary Black Hole. 

BSSN Baumgarte-Shapiro-Shibata-Nakamura 3+1 formulation of 

Einstein's equations. 

CBC Compact Binary Coalescence. 

CCATIE AEI/LSU numerical relativity code. 

EOB Effective One Body. 

EOBNR Effective-one-body waveforms calibrated to numerical data. 

ERD Effective Ringdown. 

GH Generalized harmonic formulation of Einstein's equations. 

Halmdol Numerical-relativity code developed at NASA-Goddard. 

HHT Hilbert-Huang Transform. 

IMF Intrinsic Mode Functions. 

IMR Inspiral-Merger-Ringdown. 

ISCO Innermost Stable Circular Orbit. 

LI, HI, H2, VI LIGO Livingston 4 km, Hanford 4 km, Hanford 2 km 

and Virgo 3 km gravitational-wave detectors. 

LAL LSC Algorithm Library. 

Lean Numerical-relativity code developed by Ulrich Spcrhakc. 

LazEv Brownsvillc/RIT numerical relativity code. 

LSC LIGO Scientific Collaboration. 

MayaKranc Numerical-relativity code developed at Penn State 

using the Kranc code-generation package developed at AEI, 

Southampton and Penn State. 

MCMC Markov-Chain Monte-Carlo. 

NINJA Numerical INJection Analysis. 

NR Numerical Relativity. 

PDF Probability-density Function. 

PN Post-Newtonian. 

PU Numerical-relativity code developed by Frans Pretorius. 

S5 Fifth LIGO science run. 

SNR Signal-to-noise ratio. 

SPA Stationary Phase Approximation. 

SpEC Spectral Einstein Code developed at Caltech and Cornell. 

TT Transverse- Traceless. 

UIUC University of Illinois at Urbana-Champagn numerical 

relativity code. 

VSR1 Virgo science run 1. 

WRD Weighted ringdown. 
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